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Summary 


This report presents new theoretical results which lead to new algorithms for the 
computation of fuel-optimal multiple-bum orbit transfers of low and medium thrust. 
Theoretical results introduced herein show how to add bums to an optimal trajectory and 
show that the traditional set of necessary conditions may be replaced with a much simpler 
set of equations. Numerical results are presented to demonstrate the utility of the 
theoretical results and the new algorithms. 

Two indirect methods from the literature are shown to be effective for the optimal 
orbit transfer problem with relatively small numbers of bums. These methods are the 
Minimizing Boundary Condition Method (MBCM) and BOUNDSCO. Both of these 

methods make use of the first-order necessary conditions exactly as derived by optimal 
control theory. 

Perturbations due to Earth’s oblateness and atmospheric drag are considered. 
These perturbations are of greatest interest for transfers that take place between low Earth 
orbit altitudes and geosynchronous orbit altitudes. Example extremal solutions including 
these effects and computed by the aforementioned methods are presented. 

It is a commonly accepted notion in the field of optimal orbit transfer that the 
more bums an optimal transfer executes, the lower the cost. Unfortunately, many 
numerical methods are not robust enough to simply "jump" from an N-bum solution to an 

N+ 1 burn solution - A ncw algorithm is presented which greatly eases this process. The 
method is just as easily implemented in the framework of MBCM as BOUNDSCO, any 
indirect method, or a hybrid method. 



Using this algorithm and the indirect methods mentioned above, the phenomena 
of multiple solutions is demonstrated for the optimal orbit transfer problem. A simple 
empirical guideline is proposed which chooses between two or more multiple solutions 

when using this algorithm. It is not claimed that the algorithm will obtain the globally 
optimal solution. 


Intuitively, one might want to think of an optimal multiple-bum transfer not as 
one large trajectory, but as a sequence of optimal one-bum transfers between transfer 
orbits that are optimally chosen. For ideal gravity, a strong relationship is shown to exist 
between these two problems. Based on this relationship, two new numerical methods are 
presented which iteratively compute optimal orbit transfers. The first method, named the 
Patched Method, appears to be very robust yet sluggish in convergence. The second 
method, named the Modified Patched Method (MPM) seems somewhat less robust but 
much faster in convergence. For optimal orbit transfers in ideal gravity with large 

numbers of bums, MPM seems to be superior to the other methods investigated in this 
report. 

Finally, an investigation is made into a suboptimal multiple-bum guidance 
scheme. This scheme is, in fact, seen to have somewhat less than desirable terminal 
error. This terminal error is improved through a time-to-go indexing scheme. Future 
directions for multiple-burn guidance are suggested. 

The FORTRAN code developed for this study has been collected together in a 
package named ORBPACK. ORB PACK and a user manual are provided. The manual is 
included as an appendix to this report. 
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Section I 

The Orbit Transfer Problem 

LL Introduction 

The most popular motor today for performing orbit transfers is of high thrust and 
usually a solid, sometimes a liquid rocket motor. These typically have a specific impulse, 
or I sp , in the lower hundreds of seconds (250s-450s) and thrust in the thousands of 
Newtons 1 and up. In this range, they can be considered impulsive 2 , applying changes in 
velocity on a time scale much shorter than the orbit period. For many years the study of 
optimal orbit transfer has focused on these impulsive motors. 

With the hopes of lower fuel consumption due to an l sp typically in the thousands 
of seconds, electric propulsion has recently grown in popularity and many studies have 
been performed to develop the motors; a major satellite manufacturer is already designing 
satellites which use a Xenon Ion Propulsion System (XIPS) 3 . The thrust produced by 
these motors is in the tens to thousandths of Newtons; for example, XIPS produces 18 
thousandths of a Newton with an I sp just under 3,000 sec. Obviously, orbit transfer 
maneuvers with such electric propulsion will take more time and practical transfers can 
no longer be modeled as impulsive. Since it is necessary to specify the maneuver with 
continuous functions as opposed to discrete impulsive events, the optimal transfer 
problem has been too complicated for exact analytical solutions. 


1 Hertz, J. R.., and Arson, W. J., Space Mission Analysis and Design, Kluwer Academic 
Publishers, Boston, 1991. 

2 Robbii^H^L/‘An Analytical Study of the Impulsive Approximation,” AIAA Journal, 

3 Christensen, R. A., ed., “Space Propulsion’s Latest Thrust,” Vectors, Vol 37 No 1 
1995, Hughes Electronics, Los Angeles. 
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Numerical methods for the computation of optimal orbit transfers have been 
widely studied. These numerical methods fall into three categories: direct, indirect, and 
hybrid methods. Direct methods parameterize the thrust program and then attempt to 
optimize these parameters while satisfying boundary conditions. Indirect methods 
employ the mathematics of optimal control to formulate a Two-Point Boundary Value 
Problem (TPBVP) which can then be approached with a variety of numerical methods. 
Hybrid methods are a combination of the two. These methods are often formed by 
simply removing difficult conditions from the TPBVP and optimizing some equivalent 
cost function over the released parameters. 

The main objective of this research was the computation of fuel-optimal low and 
medium thrust orbit transfers. Here, medium thrust was taken as 1 > T/W c 2 0.01 and 
low-thrust as 0.01 > T/W 0 Z 0.001. This particular definition has been made because it is 
the initial acceleration which the rocket motor produces compared with the gravitational 
acceleration at that point that determines how easily changes in the initial orbit will be 
made. In contrast, comparing the initial rocket motor acceleration with the weight of the 
spacecraft as it would measure on the planet’s surface does not directly indicate the 
motor’s ability to move the spacecraft away from a very high orbit. 

Of the utmost interest was the ability to compute highly efficient transfers for the 
ideal case. This will provide mission planners with the ability to compute a "best" 
transfer which can be used to judge more practical schemes. However, the ideal case 
does not quite represent reality; the ability to handle orbit perturbations is desirable as 
this would produce more realistic "best" transfers. For trajectories that spend much time 
near or beyond geosynchronous orbit, the dominant orbit perturbations will result from 
either Earth oblateness effects or atmospheric drag. 3 
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Software using multiple-point shooting and modified-shooting techniques were 
used and produced many solutions. Using these, some characteristics of the solution have 
been observed and studied. Identification of these characteristics has resulted in the 
development of a new method for improving optimal orbit transfers. The method 
introduces additional bums to optimal ideal-gravity orbit transfers using an under- 
exploited property of the switching function. A set of improved transfers were 
constructed and these uncovered new properties of optimal transfers. 

Furthermore, two new methods have been developed. The first is a new hybrid 
approach called the Patched Method. This method combines the robustness of a direct 
approach and the greater convergence speed of the multiple-shooting approach in a 
configuration that can handle transfers with large numbers of bums. However, the 
Patched Method pays for its robustness with speed. 

The second new method is the Modified Patched Method (MPM). MPM trades 
back some of the sluggishness of the Patched Method for a small loss in robustness. This 
trade-off is accomplished by making use of properties specific to the orbit transfer 
problem. Some of these properties appear to be new, developed here for the first time. 
Overall, MPM seems to be superior to any of the other methods applied in this report. 

The other objective of this research was the examination of a capable guidance 
algorithm for multiple-burn orbit transfer. Work on this has produced a one-bum 
guidance algorithm using neighboring optimal feedback control. This guidance algorithm 
could be used on a bum-by-bum basis to produce a sub-optimal trajectory. 

L2, Orbit Transfer MniMtnp 

The spacecraft is represented by a point mass and assumed to be a thrusting craft 
acted upon by the aerodynamic drag and oblate-body gravity forces of a central body. 
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The central body, or planet, is also represented as a point mass positioned at its own 
center of gravity. Furthermore, the problem is restricted to crafts of mass much smaller 
than that of the central body; therefore, the planet is assumed fixed in inertial space. This 
inertial space is described with a rectangular Cartesian inertial reference frame (OXYZ). 
The centra] body is fixed at the center O of this frame and the z-axis is perpendicular to 
that body’s equator. All motion within this frame agreeing with the above assumptions 
must satisfy Newton’s Second Law: 


= ^ ( 1 . 1 ) 
dt 

where m is the spacecraft s mass, v is its velocity with respect to the reference frame, and 
XF represents the sum of forces on the craft. 

In this case, gravity, drag, and thrust make up the total force acting on the craft. 
This gives 


m\ = Te T - T drag -F gravity (1.2) 

in which the thrust is some time-varying function T(r) independent of a time-varying 
direction This is most clearly derived by considering a momentum balance of the 

spacecraft as it expells mass to produce thrust; absorbing the dm/dt term into the thrust 
term produces Equation (1.2). 

The thrust direction is expressed as the unit vector ej<r). For a three-dimensional 
thrust vector the control requires a magnitude and three components or two angles. For 
two dimensional problems, the one magnitude and only two independent control 
components or one angle is required. 

It is assumed that the fuel consumption of the motor is represented by 


4 



T 

m 

gj, 


(1.3) 


where go is Earth's gravitational acceleration at sea level and /„ is the motor's specific 
impulse. 

It is assumed that the atmosphere surrounding the central body can be described 

by an exponential model as in the standard atmosphere" resulting in the following 
aerodynamic drag force: 


where fi is a constant from the atmosphere model describing air density variation in the 
prescribed altitude region. p„ is the atmosphere density at the altitude r„, 5 i s the cross- 
sectional area of the craft, C D is the craft's drag coefficient, v is the magnitude of the 
velocity v, and r is the magnitude of the position vector r. 

The gravitational potential energy to the second harmonic is 5 


_ um 




where R, is the equatorial radius of the central body. Bis the latitude angle of the current 
posrtion from the equator, p is the gravitational constant for the central body, and J 2 is a 
constant describing the mass distribuuon of the central body; for Earth y 2 -108161*10-« 
There are additional mass distribution terms, but the series is truncated here. 6 is 

"AndersonJ. D„ Fundamentals of Aerodynamics, New York: McGraw-Hill Book Co.. 
>Space Techndog X Labomrories, Flight Performance Handbook for Orbital Operations. 
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described with Cartesian coordinates by z=r cos(6). This gravitational potential exerts 
the following force on the spacecraft: 


-^-^1 + 


R ~ 


dr 

w’here N =diag {1,1,3} and I is the identity matrix. 


I r 


The equations of motion for the spacecraft are 

x(0=f(x(/)T(r),e r (r)) 

where 

x(r) = [r T (t) v T (r) m(/)] T 
and 


(1.6) 


(1.7) 


(1.8) 


f(x(r),r(0.e r (/» 


V 



r 


2 m 


e-H'-’.fSCoW 


(1.9a) 

(1.9b) 

(1.9c) 


The thrust magnitude has both an upper and a lower bound. The upper bound is 
called T max , the lower bound is zero. Therefore, the following inequality constraint must 
be satisfied for all time t e [0,ry] : 


0 < T < T, 


max 


( 1 . 10 ) 
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For the purposes of this study a simple atmosphere model was chosen. The model 
was not intended to accurately represent the Earth’s atmosphere, or any other planet for 
that matter. It is implemented only for the purpose of demonstrating the methods used 
herein and to allow examination of its effects on the optimal transfer. 

The model was defined from a reference altitude of 450 km above the planet’s 
equator. The entire atmosphere region was assumed isothermal with a temperature of 
1.000K. The density at the definition altitude was defined to be 1.184x10-12 kg/m3. The 
definition point for this model was taken from the 1976 U.S. Standard Atmosphere 6 . 
Also, it was assumed that C D - 2, a common approximation for spacecraft 7 , and the cross 
sectional area of the satellite was arbitrarily chosen to be 4 n m 2 . 

For problems in which the ideal gravity assumption is acceptable, coasting 
trajectories are well understood and can be analytically represented. Therefore, it is 
simplest to optimize the exit, or “thrust on,” point on the initial orbit and the entry, or 
thrust off, ’ point on the final orbit. A real spacecraft implementing the orbit transfer 
could simply wait in the initial orbit until arrival at the initial orbit exit point, indicating 
that the maneuver should begin. 

Hence, the boundary conditions must determine all orbital elements except 
position on orbit, and are written as 


v(*(0)-«. ( 1.1 la) 

'•'(*('/)) = «/ (1.11b) 

where the function y determines these orbital elements for the state in question and a 0 
and a f are vectors containing the desired values at the initial and final points, 
respectively. Such a determination could be accomplished several different ways. 

6 Umted States. COESA. U.S. Standard Atmosphere, 1976, Washington: GPO, 1976. 
King-Hele, D. Theory of Satellite Orbits in an Atmosphere , London, Butterwoiihs, 1964. 
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However, using the angular momentum and eccentricity vectors is perhaps the simplest. 8 
For planar transfers, all motion can be placed in X-Y plane and the components of the 
function are 

= h = xv - yu 

V'a = = [( v2 " M A) x ~ ( rT '>] (1.12) 

^3 = He y = [(v 2 -p/r) y - (r T v)v] 

Where h is the angular momentum, e x is the X-component of the eccentricity vector, and 
€y is the Y-component of the eccentricity vector. 

In the three-dimensional case, these vectors will compose six components. Since 
the angular momentum and eccentricity vectors are always peipendicular, one of these 
components will be redundant and thus removable. There is one restriction on which 
component is removed; it can be seen clearly by considering the property that the vectors 
are always orthogonal, expressed as 


V, + h ) e ) + h t e t = 0 (1.13) 

A component of one of the two vectors can be removed if it can be computed using 
Equation (1.13). In other words, since Eq. (1.13) always holds, knowledge of the 
removed component is implied and it is unnecessary to explicitly compute it. Another 
way to state this is to say that the six components are linearly dependent. Therefore, if 
for the orbit transfer problem in question, h 2 * 0 on a terminal orbit, then the \\f function 
components can be written as 


8 Kaplan M. R Modern Spacecraft Dynamics and Control , New York John Wiley & 
Sons, 1976. 
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V'l = A, = yw - ZV 

(1.14a) 

Vi = h y = zu - xw 

(1.14b) 

V'j ~ h, = xv - yu 

(1.14c) 

V 4 = = [(v 2 - /i/r)x - (r T v)u] 

(1.14d) 

Vs = P« y = [(v 2 - At/r) y - (r T v)v] 

(1.14c) 


where x,y, and z are the components of r in OXYZ and u,v, and w are the components of 
v in OXYZ. 


If the initial or final portion of a transfer traverses altitudes where ideal gravity is 
not a valid assumption, then the boundary conditions likely need to be reformulated. For 
example, a trajectory that begins at a very low Earth-altitude cannot truly coast with zero 
cost because energy would be lost due to atmospheric drag. For such a transfer, it would 
be more realistic to fix the initial point. Likewise, some missions may be more interested 
in delivering the spacecraft to a specific point in space, in which case the final condition 
should be a rendezvous condition. 

Anticipating numerical applications, note that the problem can be 
nondimensionalized. This aided by making all states roughly the same order. In the 
presentation of example solutions, the hat f) notation will be dropped and solutions are 

assumed nondimensionalized unless stated otherwise. The non-dimensionalizations 
follow: 

r = r/r* rhsm/m* (1.15a-b) 

t = (1.15c) 

and they require the following: 

vs v/V/i/r* '/ S '//V r *7i“ (1.15d-e) 
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(1.15f-g) 
(1.15h-i) 
(1.1 5j-k) 


f .= r »/ r * pepr* 

{pJc B ) b p.SC,{r*/m*) (*./„) 5 

f e (r/rn^j/Jp/r* 2 ) R, *R t /r* 

The choices of r* and m * are completely arbitrary. However, it needs to be said that 
after a problem is solved by these nondimensionalizations rescaling must be exercised 
with caution; rescaling changes the atmosphere model and changes the equatorial radius 
used for the oblateness terms. For example, a given transfer with nondimensionalized 
parameters must specify the value for R t if oblateness effects were considered. If, after 

rescaling, one intends this transfer to represent a maneuver about Earth then r* must be 
such that R, is the radius of Earth by Equation (1.15k). Similar arguments may be made 

concerning the nondimensionalized parameters for atmospheric drag effects. 

Substitution of Eqs. (1.15a-k) into Eqs. (1.9a-c) shows that the nondimensional 
dynamic equations are equivalent to Eqs. (1.9a-c) with /t=l (the value of J 2 , however, has 
no dimensions and is not changed). In Eq. (1.9a), choosing the scalings for r and t, 
shows that the only consistent scaling for v is Eq. (1.15d). Then, in Eq. (1.9b) it is clear 
that Eqs. (1.15a-h) and (1.15j-k) are required for consistency. Substitution into Eq. (1.9b) 
also shows that the factor /i appears on both sides of Eq. (1.9b), in the numerator of even- 
term; therefore, it may be dropped from both sides. Finally, substitution into Eq. (1.9c) 
reveals that Eq. (1.15i) is required for consistent scaling. 
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Section ii 


Computation of Optimal Orbit 
Transfers 

II.l. Literature Review 

One of the earliest and most notable applications of the calculus of variations to 
the orbit transfer problem was by Lawden 9 . His work established the now widely-used 
pointer vector theory. Lawden also derived many useful analytical results including an 
analytical solution for the Lagrange multipliers over coast arcs in ideal gravity 10 ; his 
expression is easily configured to trajectories where the transfer time is unconstrained. 
He went on to conclude that for the case of escape from a circular orbit, tangential 
thrusting would be nearly optimal 11 ; however, he noted that this thrust program may not 
fare so well in other cases. Lawden studied the possibility that, in addition to arcs of 
maximum thrust and null thrust, arcs of intermediate-thrust may exist in an optimal 
transfer 12 . He later wrote a general review of rocket trajectory optimization 13 and stated 
that issue of the existence of intermediate-thrust arcs was still unresolved. 

After Lawden's formulation was published, many other researchers produced 
solutions to the Lagrange multipliers over coast arcs in ideal gravity. A set of 


Lawden, D.F Optimal Trajectories for Space Navigation, London, Butterworths, 1963. 
lu Lawden, D. F„ Fundamentals of Space Navigation,” Journal of the British 
Interplanetary Society , Vol. 13, No. 2, 1954, pp. 87-101, 1954 

Lawd f n i?c 0 R ’ “9Pfl n Sf scape from a Circular Orbit,” AstronauticaActa, Vol. 4, No. 
1^5, PP* 

12 Lawden, D. F., “Optimal Intermediate-Thrust Arcs in a Gravitational Field ” 
Astronautica Acta , Vol. 8, No. 2, pp. 106-123. 

13 Lawden, D. F„ “Rocket Trajectory Optimization: 1950-1963 f Journal of Guidance 
Control, and Dynamics, Vol. 14, No. 4, 1991, pp. 705-711. 
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expressions derived by Danby 14 appear to be the earliest such work. This was actually 
for the equivalent problem of determining the matrizant. At almost the same time, Pines 
published work which derived constants of integration 15 , some which apply during any 
port of the trajectory, even intermediate-thrust arcs, and some in restricted cases. Later, 
both Eckenwiler 16 and Hempel 17 produced formulations valid in a two-dimensional 
system. Lion and Handelsman 18 derived equations for a three-dimensional system. 
Glandorf 19 produced a very useful form for the Lagrange multiplier’s that used the 
current radius, velocity, and angular momentum vectors as reference directions. Vinh 20 
developed equations which reduced the solution of the Lagrange multipliers for any 
central force field to a problem of simple quadratures. 

These analytical results have all proved useful in many studies of optimal orbit 
transfers. However, to date no closed-form expressions have been obtained for optimal 
orbit transfers, including the fuel-optimal thrust-limited case considered in this report. 
Therefore, numerical methods are used to produce exact solutions for this problem which 
has challenged the most sophisticated algorithms. These methods are traditionally 
divided into three types: indirect, direct, and hybrid. 


)4 Danby, J. M. A, “The Matrizant of Keplerian Motion,” A/A4 Journal, Vol 2 No 1 
1964, pp. 16-19. 

15 Pines, S., “Constants of the Motion for Optimum Thrust Trajectories in a Central Force 
Field,” AIAA Journal , Vol. 2, No. 11, 1964, pp. 2010-2014. 

16 Eckenwiler, M. W., “Closed-Form Lagrangian Multipliers for Coast Periods of 

Optimum Trajectories,” AIAA Journal, Vol.3, No. 6, June 1965, pp. 1 149-1 151. 

17 Hempel, P. R., “Representation of the Lagrangian Multipliers for Coast Periods of 
Optimum Trajectories,” AIAA Journal, Vol. 4, No. 4, June 1966, pp. 720-730. 

18 Lion, P. M., and Handelsman, M., “Primer Vector on Fixed-Time Impulsive 
Trajectories,” AIAA Journal, Vol. 6, No. 1, 1968, pp. 127-132. 

19 Glandorf, D. R., “Lagrange Multipliers and the State Transition Matrix for Coasting 
Arcs,” AIAA Journal, Vol. 7, No. 2, 1969, pp. 363-365. 

20 Vinh, N. X., “Integration of the Primer Vector in a Central Force Field,” Journal of 
Optimization Theory and Applications, Vol. 9, No. 1, 1972, pp. 51-58. 
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n.1.1. Indirect Methods 

Indirect methods conven the optimization problem into a TPBVP though optimal 
control theory. The most popular indirect methods by far seem to be the shooting and 
multiple-point shooting methods. 

Among the studies using indirect methods, the work by Brown, Hanold, and 
Johnson 2 ! produced an indirect method named OPGUID/SWITCH which handles 
rendezvous trajectories or free entry/exit points and free final time using a modified set of 
boundary conditions. Results with OPGUID/SWITCH were presented for medium thrust 
levels and two to three bums. 


Another indirect method, developed by McAdoo, Jezewski, and Dawkins 22 and 
dubbed OPBURN, was actually a combination of two approaches. The first 
approximated ideal gravity using a model for gravitational accelerations linearly varying 
with altitude. This assumption results in a linear steering law and was used to simplify 
low-accuracy calculation of the transfer. The data from this approach were used as the 
starting iterate of another, more accurate code. Results with this method were presented 
for medium thrust acceleration levels and two to three bums. 

Edelbaum, Sackett, and Malchow 23 produced computer code to solve minimum 
time transfers (one bum) using equinoctial orbital elements as state variables. Constraints 
on exposure to solar radiation were considered. This method relied heavily upon the 


21 Brown, K. R., Harold, E. F and Johnson, G. W., “Rapid Optimization of Multiple- 
Bum Rocket Rights, NASA CR-I430, Sept., 1969. P 

McA fWmInuf Jezewskl ; D : J >. Dawkins, G. S., “Development of a Method for 

yP •^™ aneuver Ana] y sis of Complex Space Missions,” NASA TN D-7882 
April, 1975. * 

^Edelbaum T N. Sackett, L. L and Mdchow H. L„ “Optimal Low Thrust Geocentric 
Transfer AIAA Paper 73-1074, Proceedings of the AIAA 10th Electric 
Propulsion Conference , Lake Tahoe, Nevada, November 1973. 
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method of averaging and was named SECKSPOT. Horsewood, Suskin, and Pines 24 
modified SECKSPOT to produce a code for the optimization of multiple-bum rendezvous 
orbit transfers with plane changes between circular orbits with low-thrust in an ideal 
^avity field. The transfer times for these trajectories were fixed. 

A study by Redding 25 handled point-to-point, or rendezvous, low-thrust transfers 
with plane changes. The method presented in the study includes the reduced set of 
boundary conditions established earlier by Brown, et. al. 21 It was limited to transfers to 
geosynchronous orbits in an ideal gravity field and no results are discussed for elliptical 
terminal orbits. Solutions with low-thrust were obtained for transfers with two to six 
bums. 

n.1.2. Direct Methods 

The most common technique for direct methods is to discretize the control and 
possibly the state, then optimize the performance index by varying the control and state at 
each node of the independent variable. This optimization is usually subject to some 
constraints. In orbit transfer optimization, it obviously makes sense to use any helpful 
results from the application of optimal control theory. Almost universally, direct 
methods for orbit transfer optimization make use of a bang-bang assumption which 
eliminates the possibility of intermediate-thrust arcs. The control is then taken as a 
combination of switching times and directions. 

The Direct Collocation with Nonlinear Programming (DCNLP) technique makes 
use of polynomial approximation to both perform integration and approximate the control 


24 Horsewood, J.L., Suskin, M.A., and Pines, S., “Moon Trajectory Computational 

Capability Development,” NASA Lewis TR-90-51 , Cleveland, Ohio 44135, July 


^Redding, D.C., Optimal Low-Thrust Transfers to Geosynchronous Orbit 
Lewis SUDAAR 539, Cleveland, Ohio 44135, Sept. 1983. 


,"NASA 
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at nodes. Dickmanns and Wells 26 made a significant contribution using a DCNLP 
method based on piece-wise Hermite polynomial approximations for the state and 
Lagrange multipliers. More recently, Hargraves and Pans 22 used this technique in their 
OTIS (Optimal Trajectories by Implicit Simulation) program. The Direct Transcription 
and Nonlinear Programming (DTNLP) technique is very similar to DCNLP, with 
transcription replacing collocation for implicit integration. 

Using DCNLP once then DTNLP later, Enright and Conway 2 **. 2 ? examined 
circular, point-to-point planar transfers with ideal gravity. The methods demonstrated in 
these studies were shown effective for two- and three-bum trajectories. In using the 
DTNLP method, a technique was developed for calculating the Lagrange multipliers so 
that Pontryagin’s Minimum Principle could be checked. In some cases, it was found that 
this principle had been violated. 


Vulpetti and Montreali 30 used nonlinear programming to optimize transfers 
between circular orbits with inclinations. They did include oblateness and drag in their 
gravity model; their thrust acceleration level was about 0.0019*. Example transfers 
included from two to four burns. Pourtakdoust and Jalali 31 used DTNLP for three- 


Dickmanns, F.D and Well, K.H., Approximate Solution of Optimal Control Problems 
Using Third Order Hermite Functions,” 1FIP-TC7, VI Technical Conference on 
Optimization Techniques , Novosibirsh Springer, 1974. 

Hargraves, C.R., Paris, S.W., Vlases, W.G., ‘‘OTIS Past, Present, and Future ” 

ffifton^Head °S C *1 992^ conference of Guidance, Navigation, and Control , 

28 Enright, P.J. and Conway, B.A., “Optimal Finite-Thrust Spacecraft Trajectories Using 
Collocation and Nonlinear Programming,” Journal of Guidance, Control, and 
Dynamics, Vol. 14, No. 5, 1991, pp. 981-985. 

29 Enright, P.J. and Conway, B.A., “Discrete Approximations to Optimal Trajectories 
Using Direct Transcription and Nonlinear Programming,” Journal of Guidance, 
Control, and Dynamics, Vol. 15, No. 4, 1992, pp. 994-1002. 

Vulpetti, G and Montereali, R.M., “High-Thrust and Low-Thrust Two-Stage LEO- 
LEO Transfer Acta Astronautica, Vol. 15, No. 12, 1987, pp. 973-979 (84-354) 

SlPomakdoust. S.H. and Jalali, M.A., “Optimal Three-Dimensional Oibital Transfer 

E " gineering SysKms Design an d Analysis, 
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dimensional two-bum transfers with a medium thrust level. All these studies mentioned 
above either used fixed final time, fixed entry/exit positions in orbits, or both. 

Another direct method that is gaining in popularity makes use of a technique 
called differential inclusion.32 Coverstone-Carroll, V. and Williams, S.N.33 use d 
differential inclusion concepts in a direct optimization scheme that produced one- and 
two-bum planar interplanetary rendezvous trajectories. The title of the study states that 
these trajectories are for low-thrust, but the thrust levels fit in the medium thrust range 
defined for this report. 

n.1.3. Hybrid Methods 

Methods are called hybrid if they don’t fit neatly into either of the above 
categories. Typically, hybrid methods for the orbit transfer problem involve some use of 
the Lagrange multipliers and the Euler-Lagrange equations but also use direct 
optimization to determine other parameters of the trajectory. 

Zondervan, Wood, and Caughey 34 used a hybrid method to study three-bum 
transfers with plane changes in ideal gravity and for thrust levels in the medium and low- 
thrust range. Their approach was to take the indirect setup and release the switching 
function constraint. The switching points were then optimized directly. 


32 Kisielewicz., M., Differential Inclusions and Optimal Control , Kluwer Academic 
Publishers, Boston, 1991. 

33 Coverstone - Can-oil , V. and Williams, S.N., “Optimal Low Thrust Trajectories Using 
Differential Inclusion Concepts,” Proceedings of the AAS Rocky Mountain 
Guidance Conference , Colorado, 1994. 

^Zondervan, K.P. Wood, L.J., and Caughey, T.K., “Optimal Low-Thrust, Three-Bum 
Orbit Transfers with Large Plane Changes,” Journal of the Astronautical 
Sciences , Vol. 32, No. 3, 1984, pp. 407-427. 
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Ilgen” used a hybrid scheme called HYTOP lo compute low-thmst transfers for 
an Orbit Transfer Vehicle (OTV) study. The HYTOP algorithm uses the fact from 
optimal control theory that the pointer vector function is continuous for the duration of 
the transfer. The pointer vector function, and only this function, is discretized into piece- 
wise linear functions. The state was represented by equinoctial orbital elements. The 

final mass was then optimized over the choice of the pointer vector function parameters 
subject to the TPB VP constraints. 

Each h)brid method is unique, these two are by no means representative of all that 

have been attempted. To date, there does no, appear to be any standard hybrid 
methodology. 



The following subsections describe work in this research effort using indirect 
methods and homotopy to compute solutions. Modified forms of both shooting and 
multiple-point shooting were found capable of computing medium thrust transfers with 
small numbers of bums and some low-, hms, transfers. In this domain, a new method for 
increasing the number of bums in a transfer was developed and is based a new property 
of the switching function. This new method was used to demonstrate that optimal orbit 
transfers may have multiple solutions. Also, when using this method there is a rule-of- 
thumb that may help compute the more efficient of the multiple solutions, thus, avoiding 
the need to compute all possible transfers and comparing the cos, directly. However, 
there is no guarantee of a global minimum. 


n.2.1. Application of Optimal Control 

For this problem the ch oice of performance index is clear 
35 Ilgen, M.R., “A Hybrid Method for Computing Optimal Low-Thrust OTV 
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J = m(t } ) 


( 2 . 1 ) 


where m{ij) represents the mass of the spacecraft including its fuel at the end of the orbit 

transfer. The intention, then, is to maximize the performance index, viz. maximize the 
mass at the end of the transfer. 


The TPBVP is constructed using the necessary conditions in the usual manner. 36 
Include the steering direction vector constraint in the Hamiltonian, which can be defined 
for the optimization problem as 


«(*(t)T(/),e, (!),>.(/)) = X T (<)f(x(l),r(;).e,(:))+ A.(e, T (l)e 7 (/)- 1) 

W ' VV + X ' T fm e ''^ r -|!^^ 
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(2.2a) 

(2.2b) 


>0 sp 


from which the Euler-Lagrange equations are obtained as ODEs governing the Lagrange 
multipliers 
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(2.3b) 


36 Bryson, A.E and Ho, Y.-C., Applied Optimal Control, New York: Hemisphere 
Publishing Corporation. 
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(2.3c) 
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dm m 2 
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The next Euler-Lagrange equation is easily derived as 


dH 

dt T 


_d 

de 


+ A,(e r T e r - 1} • •] = -XJ + IXfij = 0 
t Vm ' ) m 


(2.4) 


so that the necessary condition is satisfied if e r = X,/|?i v | and A, = (T|A Y |)/(2m); in other 
words, the thrust direction is parallel to A v , which Lawden thus referred to as the pointer 
vector. This choice is further supported by a sufficient condition; note that 


dH_ 

de T de T 


= 2A,I = -|A ¥ |I>0 

m 


(2.5) 


w-hen |A,| > 0, T> 0, and m finite. Also, note that if any one of these is violated during a 

bum, the trajectory is immediately indeterminate. The choice for the Lagrange multiplier 
has been made and does not need to be solved for. 


The switching function is derived by an application of the maximum principle. 
The thrust magnitude, which has bounds and 0, will give H its maximum value if it 
is at its maximum value when Hj > 0 and at its minimum when Hj < 0. The switching 
function is 




M.J L 


(2.6) 


and the switching law is 


"r>0, 7 = 7^ 

h t < o, r= o 


(2.7) 



If Hj were to be zero for a finite time the control would be singular. Higher-order 
derivatives of H T would then be needed to calculate T. In subsection II. 1., it was noted 
that this singular control has been investigated by many different researchers but no 
conclusions are widely accepted as to when, or if, it will be part of the optimal control. 

Many authors 21 - 34 - 25 - 37 have identified the switching law, and associated 
switching function, as a source of strong sensitivity in numerical solutions. 

To complete the TPBVP, the methods of optimal control supply a set of natural 
boundary conditions 


-it 


M'/H 



(2.8a) 

>•(>.)=- 


1 T 

— 

(2.8b) 

where G is defined as 




c ( x ( 7 <,)- x (^).v e) v / ) = w| 

('/) + v /[v(x(> / ))-0/] + v. T [v(x(<.))-o.] 

(2.9) 


and v( x ) w 2 s defined in Equations (1.12). Therefore, the natural boundary conditions 
can be expressed as 


37 Chuang, C.-H. and Goodson, T.D. “Optimal Trajectories of Low- and Medium- Thrust 
Orbit Transfers with Drag and Oblateness,” Submitted to the Journal of the 
Astronautical Sciences. 


20 



where 


( l / ) \ d\y / i \\1 T 

M*,)JT*M''))J v 

MOl r^/ / ui T 

moH* wo) J 1 


T W [2rv T ~ (r T v)l - vr T ] T 

** [_v x ] (v T v)l-vv T + p^ 7r (rr T -(r T r)l) 


and the subscript X denotes the skew symmetric matrix representation of the cross 
product. 

The last condition deals with the final time. For free transfer time the 
transversality condition must be satisfied 


w M'/M'/W'/))=-fp=o 


( 2 . 12 ) 


n.2.2. BOUNDSCO 


One method used here to solve the TPB VP is a modification of the multiple-point 
shooting method. The specific algorithms are those given by H. J. Oberle in the 
subroutine BOUNDSCO 38 , written in FORTRAN. 

The state defined for the optimal control problem differs slightly from the state 
used in BOUNDSCO. The state used in BOUNDSCO for numerical computation is 

38 0berle, H. J., “BOUNDSCO - Hinweise zur Benutzung dcs Mehrzielverfahrens fur die 
numensche Losung von Randwerproblemen mit Schaltbedingungen”, Hamburger 

Beitrage zur Angewandten Mathematik, Berichte 6 1987 S 
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z(t) = [x t (t) X t (t) x f v c T v/] T 

and includes a state denoting the transfer time, Xj , and the \ a and vy- vectors, from the 
natural boundary conditions. BOUNDSCO does not allow user-defined parameters that 
are determined in the iteration process, only functions of time; therefore, these last 
quantities must be included in the state z and specified to have zero derivatives with 
respect to time. Also, BOUNDSCO is restricted to problems with a fixed partition of the 
independent variable; therefore, the independent variable has been defined as re [0,1] 
with x-xx f . This requires that the system dynamics be properly transformed to the 
independent variable x so that 


d_ 

dr 
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'x(t)' 

X(t) 


Mr) 


s 

0 i 
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. v / . 


0 _ 


and these derivatives with respect to i are Eqs. (1.9a)-(1.9c) and (2.3a)-(2.3c). If x had N 
components, then the BOUNDSCO state, z, has 2N+2(N-2)+l components. 

BOUNDSCO addresses the switching function sensitivity problem by the explicit 
inclusion of switching points in the problem formulation. The number of switching 
points is not changed by BOUNDSCO. It iteratively drives the guessed switching points 
to be zeros of the switching function, Eq. (2.6). The user must then decide in which 
intervals to have the thrust on and in which to have thrust off. Unfortunately, with this 
scheme the switching law, Eq. (2.7), may not be satisfied and must be checked after 
BOUNDSCO claims convergence to a solution. 
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H.2.3. The Minimizing-Boundary-Condition Method 

The second method used herein is called the Minimizing-Boundary-Condition 
Method (MBCM) 39 . MBCM is a modified shooting algorithm in which the switching 
structure of the optimal control is implicit. The program checks the switching function 

and the switching law to ensure that Eqs. (2.6) and (2.7) arc satisfied at each integration 
step. 


As a modification to the simple shooting method, MBCM, expands the set of 
available solutions by removing one boundary condition while keeping the same number 
of unknowns. The choice of this boundary condition is arbitrary. With the number of 
unknowns unchanged, the solutions become a one-dimensional family. Since this gives a 
much larger set of solutions, it is much easier to solve the resulting boundary-value 
problem. Once that is accomplished, the search for the solution that incorporates the final 
boundary conditions is treated as a minimization problem. The gradient is numerically 
calculated and used to update the initial state until the last boundary condition is satisfied. 
This method is about as effective as BOUNDSCO in solving the two-point boundary- 
value problems for the solved optimal orbit transfers. 

n.2.4. Example Two-Burn Extremal 

A solution is presented in this subsection, obtained by both BOUNDSCO and 
MBCM. It is nondtmensionalized and assumes ideal gravity. The transfer is made 
between two planar, aligned orbits. The solution’s trajectory is shown in Figure 2.1. The 
transfer time has been optimized and is 19.05. The initial mass is 1.608. The initial 
semimajor axis is 3.847 and eccentricity is 0.02378. The final orbit semimajor axis is 1.5 
and eccentricity is 0.333. The product gj sp is 1.313 and the thrust level is 0.03. 


39 Chuang C.-H and Speyer, J.L., “Periodic Optimal Hypersonic SCRAMjet Cruise,’ 
Optimal Control Applications and Methods, Vol. 8, 1987, pp. 231-242. 
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Since initial altitude for the transfer is 3.905, the initial Tf\V 0 is 0.2845 and the 
transfer may be categorized as a medium thrust transfer by the definition stated earlier. 
With the initial orbit higher than the final orbit, this transfer may be viewed as an optimal 
descent transfer. However, since atmospheric drag has not been considered, it should not 
be viewed as an optimal de-orbiting transfer, where the spacecraft would be intentionally 
placed in an orbit low enough for drag to eventually destroy it. 

Two bums are used to complete the transfer. Most of the change in energy occurs 
in the longer second bum, but most of the change in angular momentum occurs in the 
first bum. 

D.2.5. Example Three-Burn Extremal Considering Perturbation Effects 

In this subsecdon, another example transfer is presented. This transfer was also 
obtained with both BOUNDSCO and MBCM. However, this is a three-bum transfer 
whose terminal orbits are not planar. The initial orbit has the same semimajor axis and 
eccentricity as the transfer from Fig. 1 except now the orbit is inclined 20°, has a right 
ascension of 13°, and an argument of perigee at 15°. The final orbit is also identical but 
inclined 1° with 0° right ascension and an argument of perigee at 0°. The thrust level and 
specific impulse are also the same. This solution includes oblateness effects but excludes 
drag effects. For the computation of oblateness effects. Earth’s value for J 2 (1082.61x10- 
6 ) was used along with 7?, =0.9696. Since this transfer is intended to be about the earth, 
r^=6578 km must be specified as it ensures the correct equatorial radius scaling. 
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Figure 2.1 Two-Bum Extremal Orbit Transfer Solution with Free Final Time. 


The trajectory is shown in Figs. 12-2.3. This is a fixed transfer time transfer with 
ty=28.75. Recall that this is a descent trajectoty; the initial orbit is higher than the final 
orbit. It is interesting to look at this transfer in terms of the normalized time. r. the 


energy, £, the angular momentum, h, the semimajor axis, a, the eccentricity, e, the right 
ascension, 12, the argument of perigee, a, and inclination, i, for certain segments and 
points on the trajectory. For the first bum Ar=0.3616, A£=-0.07760, and AA-0.6566. 
The bum ends at what would be an orbit of a=2.409. e=0.5420, 12=8.320°, <tfcl.l23°, and 
i= 1.665°. For the second bum Ar=0.1450, A£=-0.1048, and AAMU310. The second 
bum ends at what would be an orbit of <1=1.601, 4=0 3742, 12=-1.073°, a=0.3892°, and 
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i«1.202°. For the third bum Ar=0.02420, AE=-0.02101, and AA=-0.01865. The final 
mass for this transfer is 1.1656, the initial mass was 1.527. As a result of the oblatencss 

effects, this transfer has poorer performance than if it could be performed in ideal gravity, 
where it’s final mass would be 1.1659. 

If drag is considered in the trajectory, performance improves and the final mass is 
1.1663. This is consistent with a descending transfer whose final orbit is rather low. The 
altitude of perigee for the final orbit is 6578 km where drag needs to be considered; 
therefore, atmospheric drag can be used to improve performance. Obviously, with the 

consideration of atmospheric drag, this transfer could be considered as an optimal de- 
orbiting transfer. 

The loss in performance caused by the oblateness effect is expected. The terminal 
orbits have their apses aligned; since the oblateness effect causes the line of nodes to 
regress, the optimal thrust program must fight this effect to return the orientation to that 
of the initial orbit. The improvement caused by drag is also expected for this is a 
descending trajectory and drag encourages descending trajectories. 

It is interesting to note that the change in right ascension was almost exactly 
divided between the first two bums while the change in both inclination and argument of 
perigee happened almost entirely in the first bum. The change in inclination can be most 
dramatically seen in Fig. 2.3. The bum at the top of the figure is the first bum. The next 
two bums are difficult to distinguish but not very interesting from this vantage point. The 
second coasting orbit, or transfer orbit, is quite similar to the final orbit; fittingly, the 
third bum imparts the least energy of any of the bums. 
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X 

Figure 2.2 Projection into X-Y Plane of Three-Bum Transfer in Ideal Gravity 



Y 

Figure 2.3 Projection into Z- Y Plane of Three-Bum Transfer in Ideal Gravity 

This example demonstrates the ability of these methods to obtain exact solutions 
to the orbit transfer problem for nonplanar trajectories that include perturbing effects. 
BOUNDSCO typically can obtain such trajectories within the desired tolerance if given 
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the solution under ideal gravity as the initial guess. However, performance usually 
becomes unacceptable if the number of burns was increased beyond six; this is an 
empirical observation and by no means constitutes an absolute limitation of 
BOUNDSCO. There may well be certain cases in which BOUNDSCO can compute 
transfers with more than six bums quite easily; however, experience indicates that these 
cases are uncommon. 

IL3 A NfW Property of the Optimal Switc hing Function 
A very interesting property of the optimal control solution under ideal gravity is 
that the initial and final values of the switching function are equal. Even more interesting 

is that for the free transfer time problem they are both equal to zero at the initial and final 
times. 


This property may be explained with the following theorem. In the following, Cf 
denotes the set of /-dimensional vector functions that are continuous with respect to all 
arguments, vector and/or scalar, and U denotes the set of piece-wise continuous scalar 
functions with one scalar argument. 


Theorem 11.1 : Given a bang-bang optimal control problem of the form: 

«/ 

3 = JW x (0»0 + M x (0»O m ( 0]<* where £(*(O»0 € c ° and M(x(r),r)eC ,° 


and subject to the following: 

x(r) = f(x(r),r)+g(x(/),v(r),/)w(r), x(r)eC°, v(r) e C° ; 
UminSu(t)£u max ,u(l)€ll l 

ViWO)-0,V/(*(»/)) = 0, '•',(*(»,))*<?. ; 

and /y-arc free for optimization 
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and satisfying the following assumptions: 


(i) ^(x(r t ),r,) = 

(ii) [<fy,(x(r))/<?x(r)]f(x(r),r) = 0, [^ / (x(r))/<?x(r)]f(x(r),r) = 0; 

(iii) 

then, considering the usual optimal control formulation, introduction of the X(r) 
functions, and the Hamiltonian //(x(r),v(r),u(r),X(r),r) function 36 , the following 
statements are true: 

(1) The switching function, 5(x(r),X(/),r) = X(r) T g(x(r),v(r) ( r) + M(x(r),r), satisfies 
5(x(r i ),X(r i )) = 5(x(r / )A(r / )) = -L{x[t f ),t f )/ u(t f ) if and only if 

(2) If the Hamiltonian is autonomous with /, and tj fixed, then 
S(x(;JX(r,)) = s(x(i/),X( t/ )) and 

s(x(( / ),X(/ / )) = [w(*(i / ), v (r / ),u(r / )A(i / ))-L(*(/ / ))]/«(< / ). 


Proof: 


In the usual optimal control formulation, the boundary conditions at r ( and tj 
result in the familiar natural boundary conditions on the Lagrange multipliers, written as 

X(r,) = -(^/^x) T v ( 

l{tf) = (dy f /dx) T v f eR n 

which involve the constant Lagrange multiplier vectors v, € R q ' and Vy € R 92 , where R‘ 

denotes the set of /-dimensional vectors with real-valued components. Now, consider the 
dot product of X(r,-) and X(fy) with vectors callechj € R n and n 2 € R n , respectively: 
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MO T|, i = " v i T (^Vi/^)n 1 

j 

*V/) n 2 = - v / T (<?V//<?x)n2 

This shows that, at both the initial and final times, any vector in the null space of the 
relevant constraint gradient matrix is perpendicular to the corresponding Lagrange 
multiplier vector. Assumption (ii) indicates appropriate choices for rij and n 2 as 

n i = 

n i =f ( x ('/K) 

With these choices, the Hamiltonian at either terminal time may be written in the 
following form: 

//(x(0,v(0,w(0A(/),r) = [X(/) T g(x(r),v(/),r) + M(x(;),r)]«(;) + L(x(t).r) 

Statements (1) and (2) follow immediately. ■ 

The theorem is useful because it leads to a method for finding time-optimal 
extremals with additional w max arcs when Wmin= 0. Although not attempted in this work, it 
may also lead to a method for finding extremals with fewer Umax arcs. 

Applied to the orbit transfer problem with ideal gravity and free transfer time, 
condition (1) implies the switching function must be zero at the entry /exit points. A 
similar condition was successfully used in the place of Eqs. (2.10) by Brown, et. al. 21 for 
free transfer time problems in ideal gravity. In that work, however, the condition was 
used as a boundary condition in order to reduce the number of variables in the problem. 
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One may make more use of this property of equal switching function values than 
a boundary condition; it can be used to help add bums, improving the performance of 
extremal orbit transfers as shall be seen in the following subsections. 

II.3.1 Family of Extremals 

Exploitation of the property described earlier by Theorem II. 1, along with the 
favorable performance of these indirect methods allowed the study of the characteristics 
of families of solutions. Herein a family of solutions is defined as a set of solutions 
whose transfer times and numbers of bums vary but whose terminal orbits do not. The 

optimal terminal points will vary from solution to solution because they are free for 
optimization. 


Figure 2.4 displays a family of optimal transfers. Each data point in the figure 
represents an extremal orbit transfer by its total transfer time and final mass. The 
transfers are planar and the dynamics do not take drag or oblateness effects into account. 
Furthermore, then- terminal orbits are the same as for the transfer shown in Figure 2.1. 

Though this family appears quite disjointed, it is actually quite connected. These 
connections can be best seen by starting at the leftmost transfer (point (1) in Fig. 2.4) and 
tracing solutions of increasing transfer time. The solutions from point (1) to point (2) are 

the original set of two-bum solutions, obtained via homotopy and a TPBVP solver 
(BOUNDSCO and MB CM). 

At point (1) the total bum time equals the transfer time; point (1) is a one-bum 
solution. Point (2) represents a local optimum in transfer time; the Hamiltonian for point 
(2) is zero and this satisfies the transversality condition. 
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1.25 



Figure 2.4 Plot of a Family of Optimal Transfers as Final Mass versus Transfer 
Time 

As a result of Theorem II. 1, the switching function at point (2) indicates the 
existence of additional solutions. The situation is shown in Figure 2.5. Because of the 
slope of Hj and the fact that it is zero at both the initial and final times (from Theorem 
II. 1), the transfer may be extended optimally by the addition of a coast arc at the 
beginning and/or at the end of the transfer. This may seem trivial; one might observe that 
coast arcs can always be added; however, this particular situation leads to the addition of 
burns. Lawden’s solution 10 to the costates on a coast arc shows that on such an arc with 
a vanishing Hamiltonian the switching function is periodic. This means that the 
switching function, once crossing zero, must return to zero. In other words, for an n bum 
transfer like that represented by Fig. 2.5, the periodicity of the coast arc switching 
function hints at the existence of two different n+l-bum solutions and an n+2-burn 
solution; each by different additions of coast arcs. 

To optimally extend a transfer with coast arcs such that the switching function 
will again vanish, it is required that the switching function at a terminal orbit both be 
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equal to zero and have an appropriate sign for its slope: positive at the initial time and/or 
negative at the final rime. This situation can be seen in Figure 2.5 below, for the portion 
of the switching function labeled “Original Transfer.” 



T=0 


Extended Transfer*© 
x=0 ◄ 


T=1 


T=1 


Extended Transfer*® 

Figure 2.5 f Function to Create More Optimal Transfers; 

symbols ® and <D refer to points in Figure 2.4 


One may observe lhat the process does not guarantee a new bum - only a new 

coast arc. However, using numerical methods, one may discover that the bum can be 
lengthened. 


Adding the coast arc is trivial; lengthening the burn arc is not. The following 
bum-addition procedure worked well. To add a bum to an n-bum solution with optimal 
transfer time that begins and ends with a bum arc: Append a coast arc to the solution at 
the chosen dme, initial or final, making sure that states and costates are continuous. This 
is easily done by integrating forward from the fine! time or backward from the initial 
time. At both ends of the new coast arc the switching function must be zero. Use this 


33 






extended transfer as a guess for the numerical routine setup for an n+l-bum problem with 

a slightly longer transfer time. Finally, use homotopy to obtain an n+l-bum solution with 
a longer transfer time. 

For the guesses constructed in this report, the new coast arc was extended so that 
the switching became positive for a finite time. Since the thrust was set to T ^ for this 
new interval, the boundary conditions were violated and the new arc was a non-optimal 
bum because the natural boundary condition was violated. However, it was found that 
this new bum aided in the convergence of iterations. 

There are three options for creating the next transfer in the family: extend the 
transfer to right, extend it to the left, or extend it in both directions. However, because of 
numerical difficulties, this last option w-as not favored. First, consider extension to the 
right. Physically, this corresponds to adding the new bum closer to the final orbit. The 
resulting transfer is represented by point (6) in Figure 2.4. Starting with point (6), 

solutions with longer transfer times were easily found but solutions with shorter transfer 
times were not found at all. 

Now consider the second option, extension to the left. Physically, this 
corresponds to adding a bum near the initial orbit. The resulting transfer is represented 
by point (3) of branch (3-4-5) in Figure 2.5. Numerical difficulty was discovered in 
attempting to find a solution with a greater transfer time than point (3); how-ever, 
solutions with lower transfer times were found constituting branch (3-4-5). Additionally, 
note that this branch, though a branch of optimal solutions, is unfavorable when 
compared to branch (6-7) of the family. This example of multiplicity may be viewed as a 
rearrangement of the bums in the trajectory. It has not been shown analytically, but there 
is likely a connection to a similar result for non-optimal impulsive trajectories 18 . 
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By the above discussion, points (2) and (3) and (6) are, in fact, the same transfer. 
The only difference between these transfers is the addition of a coast arc, which makes no 
difference in the performance associated with the transfer. This means that the branches 
of the family are connected and these connections are as follows, with the transfer time 

increasing: (1) to (2) (which is identical to (6)) to (7); or (5) to (4) to (3) (which is 
identical to (2)). 

Figure 2.6 shows the switching function corresponding to the transfer represented 
by point (7). Compare this to Figure 2.5. The situation is repeating itself; the terminal 
switching points in Fig. 2.6 are close to zero. Clearly, one may attempt to expand this 
family of transfers from point (7). 



Figure 2.6 Switching Function of Transfer at Point 7 in Figure 2.5 
n.3.2 Multiple Solutions in the Family 

Evidence of the existence of multiple solutions was found. For a specified 
problem (including specification of the transfer time and the number of bums) there may 
exist more than one extremal transfer. Such multiple solutions are represented by any 
point on branch (3-4) and any point on branch (6-7) which have equal transfer times. 
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Conditions for multiplicity are not clear, but it is clear that solutions are not necessarily 
unique. It is also clear that one cannot say that just because the transfer time for one 
solution is longer than another, the former has a greater final mass; although this is 
typically an assumption made in the literature. 

One cannot help but wonder why the solutions of branch (6-7) are more fuel- 
conservative than those of branch (3-4). Both branches are extensions of branch (1-2), 
but the difference is where the new burn is placed. When the bum was placed near the 
initial orbit, far from the attracting body, the branch was unfavorable. When the bum was 
placed near the final orbit, close to the attracting body, the branch was favorable. A 
principle often seen in impulsive trajectories seems to carry over in some form to finite 
bum trajectories; it appears to be better to implement changes in velocity near the 
attracting body, where changes in velocity will produce large increases in the already 

large kinetic energy, as opposed to far away from the attracting body, where kinetic 
energy is lower. 

Finally, it is clear that during the burn addition process, one may control the 
placement of new bums. By tending to place new bums closer to the attracting body, 
undesirable solutions might be avoided. 

The possibility of multiple solutions was recognized by Brusch 40 for one-bum 
low-thrust transfers originating from a circular orbit. Brusch also provides some 
excellent analysis concerning this phenomenon. In this research, it was found that 
multiple solutions exist for multiple-bum low-thrust transfers originating from an 
elliptical orbit. That the phenomenon may occur for the more general case indicates that 
there are likely many cases with multiple solutions. 

40 Brusch, R.G. and Vincent, T.L., “Low-Thrust, Minimum-Fuel, Orbital Transfers ” 
Astronautica Acta, Vol. 16, pp. 65-74. 
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H.4. Conclusions 

In this section the indirect methods BOUNDSCO and MBCM have demonstrated 
the ability to solve the optimal orbit transfer problem for small numbers of bums and 
small numbers of revolutions. Particular solutions have been presented in some detail. 
These solutions demonstrate some effects of drag and oblateness on the optimal transfer. 

A new method for adding burns to time-optimal orbit transfers has been 
presented. This method is based on a newly observed property of the optimal switching 
function and a proof has been given for this property. The method has proven its 
practical utility by generating a family of solutions. 

This family of solutions is a set of fixed-time optimal transfers with identical 
terminal orbits and parameterized by transfer time. Using this family, some new 
properties of optimal orbit transfers have been seen: multiple-bum transfers are not 
necessarily unique, transfers with greater transfer time do not necessarily have greater 
final mass, and local optima do not necessarily occur at transitions between N and N+l 
bums when using homotopy to increase the transfer time. 

Addressing the inclusion of orbit perturbations, neither BOUNDSCO nor MBCM 
had difficulty obtaining solutions with atmospheric drag or oblateness terms. 


37 



Section in 


New Methods for Optimizing Orbit 
Transfers 

UI I- Introduction 

The bang-bang structure of the optima] orbit transfer solution is well-known. This 
means the optimal transfer is made up of a series of individual interior transfers between 
a sequence of orbits beginning with the specified initial orbit and ending with the desired 
final orbit. However, the fact that these transfers are, individually, optimal transfers has 
not yet been widely exploited. In this section, this notion is expressed concisely in a 
mathematical sense and shown to be quite useful for numerical methods. 

Two methods that originated with this notion are presented. First, the Patched 
Method is a hybrid method with a greatly reduced number of parameters. In fact, not 
only are the number of parameters reduced, but they are all free for optimization. 

The Patched Method also takes advantage of another simple idea: any interior 
one-burn transfer taken between two neighboring interior orbits of an A'-burn transfer 
should be easier to solve than the N-bum transfer as a whole. It then makes sense to 
consider using the orbital elements of each intermediate transfer orbit as free parameters. 
Given these parameters, the performance (final mass) is computed by solving each 
individual one-bum problem in succession. 

The Patched Method, however, pays for its robustness in speed. Therefore, it 
seems to be most useful as a way of refining and developing initial guesses for the second 
method, the Modified Patched Method (MPM). MPM is an indirect method; no variables 
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arc directly optimized. It enforces conditions necessary for the transfer to be an extremal 
solution. MPM assumes a bang-bang structure; however, as in BOUNDSCO, the Patched 
Method, and many other methods found in the literature, MPM does not enforce 
satisfaction of Pontryagin's Maximum Principle. For this problem, Pontryagin’s 
Maximum Principle supplies the switching law as Eqs (2.6) and (2.7). These methods 
only guarantee that the thrust will switch values at the zeros of the switching function, 
Eq. (2.6); they do not guarantee that the polarity will be consistent with Eq. (2.7). 
However, this turns out to be an easy condition to check after iterations converge. 

A few reasonable and common assumptions are made in both methods. It is 
assumed that the only forces on the spacecraft are ideal gravity and the thrust from the 
rocket motor. The number of arcs of maximum thrust is assumed Fixed; choosing the 
number of bums is often desirable and makes the problem easier to solve. The first and 
last arcs are assumed to be of maximum thrust; however, no generality is lost here under 
the assumption of ideal gravity. Arcs of intermediate thrust are assumed not to exist in 
the trajectory because numerical experience indicates that such arcs are rare if they exist 
at all. It is assumed that no pan of the trajectory will be rectilinear; in other words, the 
angular momentum vector never vanishes. Rectilinear trajectories are unlikely to ever be 
of interest in an orbit transfer problem and, if they are of interest, the implications of zero 
angular momentum should motivate the development of specialized software. 

m.2. The Pafehprt MPthofl 

Usually, when a hybrid method is formulated the assumption is made that the 
solution to this new problem is always a solution to the original problem. Intuitively, this 
is often easy to accept. However, it is even more reassuring to prove whatever 
equivalency exists between the original formulation and that used by the hybrid method. 
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This subsection describes the architecture of the Patched Method, explaining how 
it functions. Also, it is shown that necessary conditions from the traditional problem 
statement are, in fact, equivalent to the necessary conditions which arise from the 
optimization loop of the Patched Method. 

IIL2.L Architecture of the Method 

The architecture of the Patched Method is best described as an inner and an outer 
loop. Given a choice of orbital elements, the inner loop solves each one-bum problem in 
succession. Each one-bum transfer has its terminal points and transfer time free for 
optimization. However, the result is a suboptimal transfer; it lacks the optimal choice of 
intermediate transfer orbits. The choice of transfer orbits is made by the outer loop via 
unconstrained minimization of the complete trajectory’s fuel consumption. 

The method that has been chosen for the outer loop is the conjugate gradient 
method. Since such methods tend to have better performance if they are supplied with an 
analytical gradient, such a gradient was formulated for this case; the formulation will be 
presented in this section. The particular FORTRAN code is taken from a common 
reference 41 . 


The architecture of this method indicates a useful new paradigm for the orbit 
transfer problem. One might think of the multiple-bum transfer optimization problem as 
optimizing the fuel used by choice of the intermediate transfer orbits, expressed as 


given a 0 ,a N ,m o ,c,T; q min lJa w ,a 

" ’ V y-i ’ 


(3.1) 


41 Press, W.H., et al. Numerical Recipes: the Art of Scientific Computing, New York: 
Cambridge University Press, 1989. 
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where r /o = 0 and shall be called the transfer time function which 

computes the optimal transfer time for the orbit transfer problem defined by the initial 
orbital elements a,_„ the final orbital elements a„ the thrust level T, the initial mass m, 

and the fuel consumption rate c. In (3.1), the value for the initial mass of each bum is 
calculated knowing the transfer times for the burns before, giving an unconstrained 

minimization problem; alternatively this could have been expressed as a constraint on the 
minimization. 


In this section it will be proven that certain conditions necessary to solve (3.1) are 
equivalent to certain conditions necessary to solve the orbit transfer fuel-optimization 
problem, under cenain assumptions. It will be seen that the restrictions imposed are few 
and quite practical; however, it is not claimed that the two problems themselves are 
equivalent; this may or may not be true. Nevertheless, this paradigm has certain 
advantages. The problem expressed in (3.1) is a parameter optimization problem. If an 

expression for the transfer time function were available, this would quite likely be easier 
to solve than the TPBVP. 


Unfortunately, there are no analytical expressions or approximations for the 
transfer time function. The Patched Method must compute it numerically in the inner 
loop. The inner loop uses both Direct Collocation with Nonlinear Programming 
(DCNLP) and multiple-shooting to solve the one-bum transfer. Each time the optimal 
solution for a one-bum trajectory is required, either method may be used. For the first 
iteration, the choice is up to the user. If DCNLP is requested, the solution is found for a 
high tolerance. Once this tolerance is achieved, a multiple-shooting guess is constructed. 
Multiple-shooting is then used to reduce the error to the desired, lower, tolerance. If 
multiple-shooting was requested as the initial method and it fails, a DCNLP guess is 
constructed and DCNLP is attempted. If DCNLP is successful, then muldple-shooting is 
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used again. This structure was chosen because it was found that DCNLP was typically 
much too slow to use with each outer-loop iteration but multiple-shooting typically could 
not converge rough guesses. The failure of multiple-shooting typically occurred with the 
first iteration if the initial guess for the transfer was poor or the failure would occur if the 
outer loop took too large a step. 


ni.2.2. Using Direct Method Solutions as Guesses for Indirect Methods 

At this point, the question of converting the solution from a direct method to the 
guess for an indirect method arises (the inverse process is trivial because the solution 
obtained by an indirect method inherently contains more information). The adjoined 
performance index for they'th of N one-bum problems (/=1,...,N) is 

} , = '”-('-) + v i ,-, T [v(*,(0))-o ; .,] + v ! /[ v (x i (, / .))-o,] < 3 - 2 > 

+4> , h(0)-/» > ]+ K I (>)[f(n,(0.e, J (/))-x j (0]* 

where x/r) is the state, u/r) is the control, is the free final time (the initial time is fixed 

at 0), ccj.j and a } are the initial and final boundary parameters, yj(x) and y 2 (x) are the 
boundary constraint vector functions, m/t) is the spacecraft mass,f(x ; (r),e r; (r)) is the 

state dynamics, and rrijitj) is the performance index to be maximized. The parameter ft is 

fixed while solving each one-bum; its value is equal to initial mass constraint (m 0 ) or the 
final mass of the previous bum: 




(3.3) 


The discretized version for the same problem, divided into M nodes indexed by i 
and designed for a direct method, follows: 
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(3.4) 


7 ;=^ + %- ] T [C,(y ; ,)-a ; . 1 ] + V[C 2 (^*f)-aJ 

+a / T K . ~Pj} + 2 H,/A, (y,, ,to,, ) 

«*1 

where y 4 is the state, 0 % is the control, £](y) and ^2(y) arc the boundary constraint 
functions, AXy/.ty,) are integration constraints, my ,• is the spacecraft mass, and rn } w is 

the performance index to be maximized Assignment of P \ , in this case, is similar to Eqn. 
(3.3) as follows: 


Pj ~ m j- l.M 


(3.5) 


Since, for any \<k<M , both formulations solve the same problem with j—k, one 
can assume that J k - J k for any choice of a* and a k+] withm y _, M - /w y _j (r r(y _ 1} ) , then 


dJ, dJ, dJ 


dJ 


J —, and 


dJ ; 


dJ 


da, da/ da,,, da,./~~ dm,.,(, nH) ) ~ dm,. 1M 

best seen in the first-order changes for both performance indices: 


^ — . The implications of this are 


+v 2/-. T [v, ,(*,(<>))&, W-fo,.,] 

+v ./[v:,(\(', ))&,(', )-&»;] 

+ j[«,(»>(').«r,('U > (0)& / W 

0 

( x , (0. e Tj ( r )> *>> (0)&ry (0 - V*/ (r)] dt 
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(3.7) 



+1 ij>.i T [C.,(y >J )5y > . 1 -&»>.,] 

+a/[<5m ; , - 8p] 


i«i ,«! " 



Knowing the solutions for both optimal control problems, one can substitute for the state 
and control of the local extremals into Eqs. (3.6)-(3.7), respectively. The resulting 
equations are simply: 


Mj = - v 2 /<5cc, - Z/Sfij (3.8) 

67 , = ” T l2/-i T ^ a ;-] “ V<5a, - (3.9) 

It is now quite clear that since the gradients were surmised to be approximately equal, 
then v 2 j.]^r\ 2 j.], v^«r \ 2jy and 

A simple approach to converting a solution obtained with a direct method into an 
appropriate guess for an indirect method is now clear. One may use a direct method to 
compute T) 2 ;_i» T] 2; , and o ; , then use Eq. (2.8b) to obtain an approximation of the 

costates at the initial time. Knowing the states and the costates at the initial time, 

obtaining an approximate time history merely requires the solution of an initial value 
problem. 


ni.2.3. Gradient of the Cost Function 

For this application, the gradient of the cost is required. The cost for the entire 
transfer is 
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J«*raU ~ = ^ 2 '[ m Af( r A) -m i(®)] (3.10) 

>*J 1 J 

where the mass at the end of the y'th bum is a function of a Jt a;.;, and my.;. This is 
obviously equivalent expression to (3.1). Omitting some simple steps of calculus and 
algebra, the gradient of the cost functional Joveraiu may easily be written as 


dJ_ _ ~gj, r 

fir (*/(;♦!)) 


dm -4 { n‘*i)) ^(h) 

da, ~ T 

1 (*> ) 

da, 

dm, (i^ da, 


dJ _-ZoK 

*»w(Vr) 

^(v) 

T 

da s _ , 

‘K-i(V-i)) da ,-, 


Equations (3.11) are not yet sufficient to implement the Patched Method. 
Expressions for evaluating the terms in Eqs. (3.1 1) are required. To begin, note that my is 
the performance index of the yth bum. Referring back to Eq. (3.8), one observes that 


BJ, ,) T 


(3.12a) 

H 

it 

i 

< 

KJ 

(3.12b) 

dJ , _ 

(3.12c) 


so that Eqs. (3.1 1) can be restated as 


dJ 

da, 


T 


N - 1 


IK-U 


u* 1+1 


[ v 2 i+i T + £. +1 v 2l T ], i = 1 N-2 


dJ_ 
d&N - 1 



(3.13) 
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Which, simply, gives the gradient of the overall cost function in terms of the Lagrange 
multipliers from each respective one-bum problem. It is interesting to note that zeroing 
this gradient supplies simple relations between the Lagrange multipliers associated with 
the beginning of one bum to those associated with the termination of the previous bum. 
It is the "patching" together of optimal bums implied by these relations that inspired the 
name of the Patched Method. 


m.2.4. An Equivalent Set of Necessary Conditions 

The following results will prove useful to showing the practicality of the Patched 
Method conditions and, later, the practicality of the Modified Patched Method conditions: 

Lemma m.l: If the matrix T e R yields rank(T) = n - 1 and satisfies 

Tf = 0, f e R H while f satisfies X T f = 0, X e /?” and f T f * 0, 
then X may be expressed as X = T t v where v e R*~\ 

Proof: 

If rank{T) = n- 1, Tf = 0, and f T f *0, then f is in the null space of T and it is 
obvious that ranker'* f]) = « . This in turn implies that there exists a ve R* 1 and (3 e R 

such that 


X = [r T 



Now, X T f = 0 => v t H + pf T f = 0 => /3f T f = 0 => p = 0. 


Lemma III.2. Consider the following system of ordinary differential equations: 
(i) ~\(r)=(([) 
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(ii) 




MO 


and a matrix function r(x), if ^T(x(/)) + r J^(*(0) = 0, then the vector 
function MO = r(x(r)) v is a solution to the differential equation (ii). 

Proof: 


To show that a function is a solution to (ii), it suffices to substitute the function 
into both sides of (ii) and show that equality holds. 


L.H.S.= 
R.H.S. - 


f( r (x«) T v)=[fr W , ) )] 1 , 

- ^f(x(/))J T r(xd)) T v 


The left hand side will 


equal the right hand side if 


|r(x(,)) + r^f(x(,))=o. 


The following definitions are precursors to a theorem that will prove the 
equivalence between necessary conditions for the Patched Method, which will be 
expressed in the definition of conditions (//}, and necessary conditions derived from the 
usual application of optimal control theory, which will be expressed in the definition of 

conditions {/). The specific problem formulation for which such conditions are 
equivalent will be defined as (P). 

In what follows, Cf denotes the set of /-dimensional vector functions that are 
continuous with respect to all arguments, vector and/or scalar, and U denotes the set of 
piece-wise continuous scalar functions with one scalar argument. 
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Definition: The optima] control problem [P] is of the form: 

minimize J = >'(t / ) subject to the following constraints: 

x(r) = f(x(/))4g(y(r),v(r))u(r), x(r)eC°, v(r)eC°; 

>'(r) = cu(r), >’(r) € C,° 

0 £ U(t) £ ,u{t) € U ; 

v(x( 0 )-a. =0 , vfxy-o^O, V(x(0)eC, ; 
y(t 0 )=y 0 ; 

tj is free for optimization, t 0 is fixed 
and satisfying the following assumptions: 

(i, [|^(*«)] r (*«)= 0 ; 

(ii) and the number of arcs with is N 

(iii) g(x(r),y(r),v(r)) is not linear in v(r) 

(iv) the solution only contains arcs with u= 0 or u=u max ; 

(v) rank ~(x(r))1s n- 1 ; 

+ when *(0= f(x(r)) 

(vii) f T (x(0)f(x(0)’‘OV/€[t„t / ] 

Consider the usual optimal control formulation, introduction of the Lagrange 
multiplier functions i{t), the Hamiltonian //(x(r), >'(/), v(/) f ii(/),X(r)) function, and the 

following panition of X(r) 
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Definition: For optimal control problem (P), the conditions {/} are 


"(x(0.y(0 > v(/),n(r)i(r))=X l T (r)f(x(r)) (3.14) 

+ [*« T (0 g(y(0rV(r)) + cX,(r)]«(r) = 0 

^»(0“-[^f(x(r))j MO (3.15) 

J t K (0 * (Oj^; g(>’(0> v(r))ju(r) (3.16) 

^ (r)T [^ g(> ' (r),v{ ^] = 0 (3.17) 

^( r /) = [^( x ( r /)) v, (3.18) 

^( ,# ) B ”[^( X ( r ')) C3-19) 

*>('/) = 1 (3.20) 

K(L ) T g(>’fc, ). y(r ti )) ■ + c \ y (/„ . ) = 0, / = 1, . . . 2(N - 1) (3.21) 


These are the transversality condition, Eq. (3.14); the Euler-Lagrange differential 
equations, Eqs. (3.15)-(3.17); the natural boundary conditions, Eqs. (3.18)-(3.20); and 
that the switching function vanishes at the switching points, Eq. (3.21). It is also required 
by conditions (7) that the control u(t ) switch values across each switching point, in a 
pattern consistent with assumption (ii). 


Definition. For optimal control problem {P}, the conditions {//} arc 
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(3.22) 


^■(*(0.y(0.v(0.tt(/U < (o) = x ll .(/) T f(x(0) 

+ ^(0 T g(>’(0. v (0)+ ^,(0^(0 = o 


|^(' ) =-[£ f( * ( 0)] T x 1 ,(, ) 

(3.23) 

■£*,(')= -K(if ^e(>’(').v(l)) u(r) 

(3.24) 

“(0 = 

(3.25) 

M0 T *J~g(>’(0* v (0) =0 

(3.26) 

^.(0=-[f«0)] T v. 

(3.27) 


(3.28) 

f— « 

II 

(3.29) 

+ ^>(,4l)( r i + l)Vy! = 0 

(3.30) 

'..1 

X ('.*l ) = *('/!)+ Jf(x(0M 
1 - 

(3.31) 

fi 

' >’(0 = >’(^ 1 ) = >’fc) • 



“(0 = 0, /e^.r^] 

where Eqs. (3.22)-(3.26) are defined for/ e and the following partition is 
defined 
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All conditions in {//} arc defined for i=l..JV except Eqs. (3.30)-(3.31) which are 
only defined for »=1..JV-1. Finally, /j=/ 0 is assigned and the value for ^ is seen to 
be ip. 


Theorem 7/7.7: If and only if 

{x(0.M0.v(0,tt(0A(r)|r«[r # ,i / ]},v #f v /f r /) {^|/«i i ...2(^-l)} (3.32) 

satisfies {7} then 


{x(r),>(r),v(r),w(r)| / e [/pr A ]},{(x,(r), t [ i = l,...tf} 


(3.33) 


satisfies {//), assuming that the constraints and assumptions from {/*} are 
satisfied. 


Proof: 


It will be shown, for both the necessary and sufficient parts of the theorem, that if 
one condition holds, then a construction may be made such that the other is satisfied. 

Assume that (3.32) satisfies {/). A solution to (77) will be constructed from 
(3.32) going backwards in time. For the last u= Umax arc, where/ e [/„,/,*], define 


V A =V / 




/ € [/yv^/Tv] 


(3.34) 

(3.35) 

(3.36) 


These definitions allow Eqs. (3. 14)-(3. 18) and Eq. (3.20) to imply satisfaction of Eqs. 
(3.22M3.26), (3.28), and (3.29) fon e [/*,/,*] and i=N. Eq. (3.21) for /=2(W-1) specifies 

that the switching function is zero at the beginning of this interval, where t-t N . 
Therefore, satisfaction of Eq. (3.22) for i=N clearly implies thatX 1Ar T (/ /v )f(x(/ A ,)) = 0. 

Considering this result. Lemma III. 1 with T(x(t„)) = |i£fx(r„)) and assumptions (i). (v), 
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and (vii), implies that there exists a -v eN e such that Eq. (3.27) is satisfied for i=N. 
This completes the definitions for the final arc. 

Consider the next interval, where / e the definitions will now be 

extended into this interval. Define !-])~U(2S-3)- The conditions {/} specify that w(r)=0 
for ; in this interval. This implies that Eqs. (3.31) with i=N-l are consistent with the 
switching structure of (/). Define 

*,a(0 = M0> 

With this definition and that Eq. (3.27) is satisfied for i=N, Lemma III.2 with 
T(x(r)) = ^~( x (0) and assumption (vi) implies that the Lagrange multipliers satisfy 




dx 


The definition v /(A ,. |) = - - j - - v oS then implies that Eq. (3.30) for i=N-\ is satisfied. 
The construction for the last u=0 arc is complete. 


Define 


*A'-1 “ f»(2 \-4) 

^- i(0= rh^ (0 ' r€ [ rjv - i,r /^-»] 

Note that this definition implies satisfaction of (3.29) for i=N-l because 

This also makes satisfaction of Eq. (3.30) for j=AM imply 

satisfaction of (3.28) for i~N- 1. After establishing these constructions, the arguments for 
the previous and t/=0 arc may be repeated. With each repeat, the construction is 
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made with scaling by an even earlier value from X^t) in the following sequence 
^■,(0’ * = W....2. Such repetition may be continued until the beginning of the first bum 
is reached. At this point, the definition 

1 


implies satisfaction of (3.27) with /=1 and completes the proof of the “if’ part of the 
theorem. 



Assume that (3.33) satisfies {//). The construction of the solution to {/} will 
proceed backwards in time. Consider the last arc, where r e [w]- Define 

V/ = Vy* 

fi2(W-l) = 

= M 0. ^[r^,r A ] 

For t and i=N, this construction lets Eqs. (3.22)-(3.26) and (3.28) and (3.29) 

imply satisfaction of Eqs. (3.18) and (3.20) at the final point and Eqs. (3. 14)-(3. 17) 
during the interval. Now, it is obvious that satisfaction of Eqs. (3.14) and (3.27) with i=N 

in this interval under assumption (i) implies that Eq. (3.21) is satisfied for /=2 (jV- 1); in 
other words x, 2(A ,_,j is a switching point. This completes the construction for the last 


The definitions will now be extended into the interval With Eqs. 

(3.31), the conditions (//} specify that u(t )= 0 for t in this interval. Define 
This implies that Eqs. (3.31) are consistent with this switching structure of {/} up to and 
including this interval. Now define 
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for all i in this interval. Knowing u(t)~ 0 and that Eqs. (3.31) are satisfied in this interval, 

Lemma III. 2 with assumption (vi) and r(x(r)) = -^-(x(/)) implies satisfaction of Eq. 

erx ^ 

(3.15) in this interval. Define 

A ,(0 = ^ yS ) 

for all i in this interval. Knowing u(t)= 0, this immediately implies satisfaction of Eq. 

(3.16) in the interval. Finally, since Eq. (3.14) was satisfied in the previous interval, Eqs. 
(3.15M3.16) are satisfied continuously from t=tj to any point in the current interval, and 
since the control switched values at a switching point, then Eq. (3.14) is satisfied in this 
interval. This completes the construction for the last u= 0 arc. 

Define t s(2N ^^i s .j. Consider the interval Conditions (//] specify 

that this is a u = 1 *^ interval which, by the definitions, is consistent with the switching 
structure of {/}. Define 


MO- )^i(.v-i)(0 

*>) = M ^) V -,)(0 

in this interval. Equations (3.22) and (3.28) with i-N - 1 imply that is a switching 
point. Considering the definitions, Eq. (3.28) with i=N-l and Eq. (3.30) with i-N- 2 
obviously imply continuity of the Lagrange multipliers X,(r) across the switching point 
l J{N continuity of X y (t) across this point is immediately implied by the definition. 

Therefore, Eqs (3.15) and (3.16) are satisfied across the switching point. 



The previous arguments for the final u=u max and u=0 arcs may be repeated, 
implying satisfaction of the conditions in {/} for each interval. After repeating the 
arguments and reaching the beginning of the trajectory, the following definitions will 
have been made and are presented for the sake of clarity: 


n^yfc) ^(x(0)] [-v.,], r€ [ r /('->)4’ i = 2,...N -l 


A,(0- 


nU',)\ ' € ['/m)4 '=2,...n-i 


X(0- 


N 


m>,) 


X,(r), i= — 1 


Finally, for the first interval, one more definition is required. The definition 



forces satisfaction of Eq. (3.27) with /=1 to imply satisfaction of Eq. (3.19). ■ 

The theorem does not assure satisfaction of Pontryagin’s Minimum Principle. 
This principle requires that 


“(0 = 0 when M0 T 8(y(0.v(0) + cX,(r) >o 
u (0 = when X 1 (r) T g(y(r),v(r)) + cX > (r) < 0 

It should be noted that in the application of the Patched Method to the optimal 
orbit transfer problem, a second-order condition was taken into account. Lawden’s 
pointer vector theory is a second-order condition and is explicitly specified. Also, note 
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that this condition was determined considering the maximization problem instead of the 
equivalent minimization problem. 

To apply Theorem III. 1 to the orbit transfer optimization problem, the 
assumptions of the theorem must be satisfied. Assumptions (i), (iii), and (vii) are 
obviously satisfied. There may still be debate over assumption (iv); however, based on 
numerical experience, orbit transfers that violate (iv) are rare if they exist at all. 

Assumption (ii) is made in anticipation of the ideal gravity assumption. In such a 
case, coasting before the first bum contributes zero cost and coasting after the final bum 
contributes zero cost. It therefore makes no sense to allow such arcs as part of the 
trajectory to be calculated. If an initial and/or final coast arc is desired, it may be added 
to the computed trajectory without affecting optimality. 

Rectilinear orbits will be explicitly excluded from candidate orbit transfer 
trajectories. Such orbits intersect the center of gravitation and are, therefore, rarely of 
interest for the orbit transfer problem. With this exclusion made, assumptions (v) and 
(vi) may now be shown true for the orbit transfer optimization problem. 

It is desired that if h ~ |r x v| ^ 0, then the vector function 


( 

>1 

\ 

r x v 

< 

X 

II 

< 

z'" 

V 

= [lj*S Ojxl] 

vx(rxv) — ii-r 
Vr T r 


yields ^“ 5 - N ° te that thiS formuIation for V(*) calculates the angular 

momentum and eccentricity vectors, then removes the third component of the eccentricity- 
vector. v(x) as defined above yields 
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fW-rt 0 i 

ax “i 1 *-* Mi 



I-v.] 

[«-.] 

® 5 «l] 

[ ( V T v)l - VV T + (rT ^, , ( rr T - (fT r)l ) 

[ 2 rv T -(r T v)l- vr T ] 


where the subscript "X" denotes the skew symmetric matrix representation of the cross 

^ I - 5 is desired. The task is simplified by the 


product. The result, rank 
following simple manipulation 


Mswi.f, o i 

ax i u} U5 *>J 


^3x3 0), 3 

V X ^3x3. 


1 L (- T 0 


(- v -l [r.] 

[rv T -vr T ] 


0 , , 

T \ 3'2 r x r X 


which makes use of ihe identity a x b x = ba T - (a T b)l . This, in combination with 

I,„ 0, •>' 


rank 


[^3x5 0 Jxl ] 


'3*3 

V x ^3x3 J 


= 5 


implies 


rank 


■¥*('))] . I 

5 , rank 


l-'.l 

K] 

y 

r min \ 


- . 

[rv T - vr T ]j 

> 

\ j 


K . 

( T \3 :2 r X F X 

i( r r ) J 

A 


It is most convenient to consider, without loss of generality, the following rotation of 
vectors r and v into the X-Y plane via an orthonormal matrix W defined such that 


Wr = 


X 


u 

y 

and Wv = 

V 

o_ 


0_ 
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It is easy to show that this rotation does not affect rank 
that, after rotation, 


cfy(*(0) 

ax 


. Substitution reveals 


rank 


t" v J 


T 2 r X r X 






f 

' 0 

0 

-V 

0 

0 

y ' 

\ 

Nl ' 

> 



0 

0 

u 

0 

0 

-x 


r t t l 


= rank 


V 

-u 

0 

y 

X 

0 


[rv -vr J 




-gy 2 

gxy 

0 

0 

-A 

0 


- 

) 



gyx 

-gx 2 

0 

A 

0 

0 



\ 


0 

0 

0 

0 

0 

0 

mi 

/ 


where h=xv-yu and g - . It can be shown that 

(r r) 


del 


' 0 
0 
v 

-gy 2 


o 

o 

-u 

gxy 


L gxy -gx 


det 


' 0 
0 
v 

-gy 2 


-V 0 

u 

0 

0 -A 0 
0 0 




0 

x 0 


0 
0 

-u 

gxy 

L gxy -gx 2 


o 

-v 0 y 
u 0 -x 
0 y 0 
0 0 0 
0 A 0 


= ~gxh 3 


= g>'h 3 


<M X (0) 


so that as long as A*0, — — — has a nonzero minor of order 5. In other words, as long 


as the orbit is not rectilinear, rank 


0) N 

<?x 


= 5. 


Now, for assumption (vi) it must be shown that if the vector function f(x) is 
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V 

H _ 

r t 
(r T r) 


and y(x) is as already defined, then when ~x(r)= f(x(r)). 


(d d 


** v(xW ) + *v(»(0)£f(x«)))=o 


It is easy to show that 


s f «- 


[o] 


'^ I+3 MF rrT 

LL t r r i ( r r ) 


[i] 

[0] 


r^V)lV hC l"* n0tati ° n h3S bCCn dr ° PPCd f ° r COnvcnience ‘ Evaluating 


gives 


o s „] 

M„=-v, 

( r r ) 


M„ M„ 


M„ M 


22 . 


where 


>' = ( VT v)l-VV T + ^K-( r T r )l) 

= {-( rIy ) + 2 (" T ) + 2(vr T ) - ^ (rr T )| 




Next, the time derivative of each term in £- V (x) can be expressed as: 
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With these expressions it can easily be shown that 


This is more than just satisfaction of a simple condition that proves useful to the theorem. 
In fact, this shows that Eq. (2.12) is the solution of the ODEs for the Lagrange 
multipliers, Eqs. (2.3a-c), when the Hamiltonian vanishes and ideal gravity is assumed. 
As reviewed earlier, many previous research efforts have focused on obtaining such 
solutions, but the form found herein is different from those. 

HL2.5. Solution using the Patched Method with Eleven Burns 

The plots below represent the current capability . of the Patched Method. The 
eleven-bum solution represented by these plots has a larger number of bums than 
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obtained BOUNDSCO or MBCM 


*■ in 11115 stud y- Fcw solutions, if any, with this number 


of bums have been obtained in the literature. However, the Modified Patched Method, 

introduced in the next subsection, has produced solution with even larger numbers of 
bums. 


Also indicative of the Patched Method, the convergence tolerance for the outer 
loop was set relatively high, 10 ’, to prevent prohibitively long computation times. 

For this example, the thrust level is 0.09698, the product gj v is 0.3929, the initial 

mass is 10. The initial orbit is circular with a radius of 1; the final orbit has an 

eccentricity of 0.398 and a final semimajor axis of 1.708. With this information the value 

of T!W C for this transfer is calculated io be 0.009698, placing it in the low-thrust transfer 
range. 


Figure 3.1 is a plot the transfer orbit elements, viz. angular momentum, 

eccentricity vector x-component, and eccentricity vector y-component, versus transfer 

orbit number. The shape of the angular momentum and eccentricity x-component curves 

seem to indicate a second order polynomial fit could be used to reduce the number of 

vanabies in the problem. The eccentricity y-component is always small in this transfer; 

suggesting that it could be assumed aero or, more generally, the same parameterization 

may be used. The zeroth orbit is the fixed initial orbit and the eleventh orbit is the fixed 
final orbit. 
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Figure 3.2 shows the angular position of the initial orbit exit point and final orbit 
entry point of each versus the index enumerating which transfer orbit the bum ends at. 
The symmetry of this plot is somewhat sutprising. Even though each transfer orbit has its 
apse roughly aligned with the x-axis, each pair of angular positions are not reflected 
about the x-ax.s. The trend over time is almost exactly opposite between the two 
positions, but note that the values are not quite the negatives of each other. Also, it is 
clear that each bum of this transfer are perigee bums; each occuning around perigee. 
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Figure 3.2 


Orbit Transfer Terminal Points Indexed by Ending Orbit 


Another interesting trend is found in Fig. 3.3, showing the bum length versus the 
same index as before. The bum length decreases monotonically with each successive 
bum, but does not decrease linearly. One can, of course, observe a relationship in the 
trend of bum length and angular positions from Figure 3.2. Both plots have a shatp 
change a. the third bum which holds till the fourth bum and then returns to follow the 

trend from the first two. The irregular trend for this bum is attributed to the high 
tolerance given for the convergence criteria. 
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Figure 3.3 Transfer Time Indexed by Ending Orbit for the Eleven Bum Solution 

JTI.3. Th e Modified Patched Method fMPMl 
The Relaxed Patched Method is tailored to the orbit transfer optimization problem 
through known relations concerning the behavior of states and costates at different points 
along the trajectory. The concept central to these relations is that each bum of a multiple- 
bum orbit transfer qualifies as an optimal transfer between its own local terminal orbits. 
This method uses an algorithm similar to shooting methods. 

This method puts forth an algorithm for computing problem constraints given the 
values of the problem variables. The number of variables and constraints are equal. 
Also, the method can be used with any multi-dimensional root-finding algorithm. The 
discussion below describes the variables and computation of the constraints for a two- 
bum trajectory. 


64 



In the following description of the variables and constraints, the vector 
*■ = [*•/ *•/] is used instead of the more common [X r T \, T A„] T so that A m can be 

discussed separately. 

The arc between points #1 and #2 is assumed to be an arc of maximum thrust. 
Referring to Fig. 3.4, the variables at #1 are the initial true anomaly, 0 y ; the first bum 
length, tjj\ and, the vector of constant Lagrange multipliers for the Stan of the first bum, 
v »* The onl y constraint associated with point #1 is for v, to have unity magnitude. 



Figure 3.4 Diagram Illustrating the Layout of a Two-Bum Transfer 

Knowing the true anomaly, 6, and the rest of the orbital elements, a, state, x(t) 
may be calculated with the function x(0; a). Therefore, the Lagrange multipliers, A(r ; ), 
and the state, x(r y ), at the initial orbit exit point may be computed using 

x ( r i) = *(0i-' a o) (3.38) 
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(3.39) 


M0= |^W0) v, 

Where y(x) is a function that calculates the orbital elements a given the state x. The 
Lagrange multipliers, X(tp), and final state of the first bum, are calculated by 
numerical integration of the Euler-Lagrange and state differential equations. 

The vector variables Oj and v 2 are associated with point #2. These are used to 
evaluate the constraints at point #2 as 



(3.40) 

%')• /t-W'/.)) v * 

(3.41) 


The trajectory between points W2 and W3 is assumed to be an arc of null thrust. 
The variables 6 2 , the initial true anomaly for the second bum, and tp, the second bum 

length, are associated with point #3. With these values, the Lagrange multipliers and the 
state may be calculated, much as before, with 

x(/ 2 ) = x(0 2 ;ai) (3.42) 

M' 2 )= J^xk)) v 2 (3.43) 

Using the integration results from the first bum and Eq. (3.43), the following constraint is 
evaluated at point #3 

M r />)| = M' 2 )| (3.44) 
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The arc between points #3 and #4 is assumed to be of maximum thrust. The 
variables 0 2 , a„ and v 2 , specified at points #2 and #3 enable the calculation of the 
Lagrange multipliers, and final state, x^), in the same manner as the previous 

bum - numerically integrating from t 2 to tp with the initial conditions Eqs. (3.42) and 
(3.43). 


The two-bum trajectory ends at point #4. The constant Lagrange multiplier vector 
v 3 is associated with this point. The constraints evaluated at point #4 are 



(3.45) 



(3.46) 


These constraints complete the system. 

With the discussion of the formulation for a two-bum trajectory concluded, the 
formulation for a more general problem is clear. For an W-bum trajectory with a 0 , a A „ 
m o » 8 o' Isp specified, the variables are 


By use of w'hich, the following quantities are calculated 

x (O-*(0,. 'a,-,) ; i = 


M0= 


$<*'■» 


it 


V,; i = 


(3.47) 


(3.48) 

(3.49) 
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(3.50) 


x M B *M + fr(*(0) + -Y7v(t)Jr 

m{t) 

= + | -£f(*(i))] T X(;)A 

where v(r) = 

' |M')| 

and m (r) = m ,--I-(r-r / , + i(r jS -r / )|, /€[,„;,] 

. &o l sp V y* 1 / 

The constraints that must be then evaluated and satisfied are 




N-i 

v(x(^)) = a,; 

X (^) = [x( x fc)) 


v 1+1 ; i=\,...N 


.N-i 


(3.51) 

(3.52) 

(3.53) 

(3.54) 

(3.55) 

(3.56) 

(3.57) 


This gives a total of 2A , (M+1) variables and the same number of constraints, where M is 
the number of orbital elements. For nonplanar transfers A/=5 but for planar transfers, it is 
more efficient to rotate the coordinate system so that Af=3. 


In summary, the Modified Patched Method executes the following procedure for 
the ith bum, i=l ..TV, of an N-bum transfer. Given the current iterates 6 l , a , and v, , 
(note, however, that a 0 is not an iterate but a specified constant) calculate x(r,) and X(? ( ) 
with Eqs (3.48)-(3.49). If i=l, evaluate the scaling constraint, Eq. (3.54). Given i fl , and 
the calculated initial values x(r,), X(r,), compute x(r y7 ), Ify) with Eqs (3.50)-(3.53). 
Evaluate the bum terminal point constraints, Eqs (3.55)-(3.56). If i<N, evaluate the 
switching function constraint, Eq. (3.57), where KOi+j) is calculated with (3.49) knowing 
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When implementing MPM on a computer, the angular variable 6, should be 
replaced by the variables / J( , and the constraint / 1< J + This common substitution 

removes the periodic redundancy that may confuse a numerical method. 

Completion of the iterative process updating the variables in (3.47) to satisfy the 
conditions in Eqs. (3.54)-(3.56) allows the final condition of the Modified Patched 
Method to be checked. Briefly, this checks the switching law: 


MOl KM 
m (0 U, P 


>o.r = r_ 


Mol ^ w (o 

m (0 gj, p 


< 0 , 7 = 0 


(3.58) 


This condition is, in fact, borrowed directly from the application of Pontryagin’s 
Maximum Principle. When all conditions are satisfied, it may be claimed that an 
extremal solution has been obtained. 


The relationship between the Patched Method and MPM is primarily in the use of 
Eqs. (3.49) and (3.56), which perform basically the same function as Eqs. (3.27), (3.28), 
and (3.30) from the Patched Method. However, MPM also includes a technique 
apparently first employed by Brown, et. al. 21 which removes one Lagrange multiplier 
(A m ) and significantly affects the way the switching conditions are handled. This 
technique is present here as the use of Equation (3.57). 

m.3.1. Equivalency of MPM Conditions and Necessary Conditions 

This subsection is concerned with proving the equivalency between necessary 
conditions and the Modified Patched Method conditions. From the standpoint of showing 
mathematical equivalence, some combinations of variables and constraints in MPM are 
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unnecessary. Essentially, guessing intermediate orbital elements can be replaced by 
requiring the state to be continuous between bums. 


Definition: For optimal control problem {P}, the conditions {///) are 

l v .|>° 

MO* ^(*(0) v,;/= I...N 

^fc) T g(>^). v (^))*X(r 1+ J T g(y(r 1+ J I v(/ (+] )) ; i = 

x fc)* x (O + )[f(x(0) + g(>(0.v(0)^]d/ 

A J 1 i = 

where X(r) T ^g(>’(/),v(r)) =0 

and >W = + + 'e 

x (',*i)= *('/;) + Jf(x(r))^r 

* >’(0 = >’(^) = >'(^) I = i,...yv -1 

«(0 = o, /e[r,,r (+1 ] 

x (v) T g(>’(v)'V(v))l 


where ij=t 0 is assigned and the value for tyis seen to be 
Theorem 1112: If and only if 
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{x(^).>(r). v(/) # «(r), X(r)|/ e [/ 0 , r 7 ]}, v o , v, , r 7 , {r„ |/ = 1, . . .2(A^ - 1)} (3.66) 


satisfies (/) then 


{x(0.>'(0-v(r),«(r),?L(0|f e[r I ,f yJV ,]},{r i ,r ys (/= i,.„a^},{ V j |/ = + (3. 

satisfies (///), assuming that the constraints and assumptions from (P) 
satisfied. 
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Proof: 


Both sufficiency and necessity will be proven by assuming satisfaction of one set 
of conditions and then constructing the solution to the other. From here on, assume that 
the constraints and assumptions from (P) are satisfied. The “if’ pan will be proven after 

the “only if’ part. To prove the “only if’ pan, it will be useful to follow time in reverse 
from t=ij to the initial time, t=t 0 . 

Assume that (3.67) satisfies {///). Define a scaling factor yeR, 


MvIsMvMv)) a6S> 

Equation (3.65) ensures that the y exists as a finite real number. Define v = yv 
t( r /)= and recall that t^t^. Note that this construction makes satisfaction of 

(3.60) with WV imply satisfaction of (3.18). Now, define A, (t f ) = 1 which satisfies 
(3.20); this makes the switching function in the form of Eq. (3.21) vanish for t=t f . 

It is obvious that when assumption (i) holds, Eq. (3.18) is satisfied, and Eq. (3.21) 
vanishes for r^then Eq. (3.14) is satisfied at r*r, Now, extend the construction so that 
M0=rM'), re[r*,r A ] andEq. (3.16) is satisfied. Note that this and Eqs. (3.63) imply 
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that all Euler-La grange differential equations, Eqs. (3.15)-(3.17), are satisfied in this 
interval. Therefore, the Hamiltonian is constant in the interval and is hence equal to zero 
at t=i N . Now, with the Hamiltonian zero, assumption (i) and Eq. (3.61) with i=N implies 
that the switching function vanishes again at t=t N . Define t s2(N . } ft N . Since by (3.64) and 
(3.63), the bang-bang control, u(t), switches from u ^ to zero at t=t N , the Hamiltonian 
will be continuous across this switching point and, therefore, zero. 

Lemma ID. 2 with r(x(r))=|^(x(0) for t Eq. (3.64), and assumption 

(vi) implies satisfaction of Eqs. (3.15) and (3.17) in this interval. Extend the construction 
* * + 

so that X y (t)- A y (r A -) = in the interval, thereby satisfying Eq. (3.16). Having 

this construction, knowing that the switching function vanishes at t=t^ that u(:)=0 is 
assigned in this interval by (3.64), satisfaction of Eq. (3.62) implies that the switching 
function vanishes at In order to imply satisfaction of Eq. (3.14) at the end of this 

interval, it must be recognized that again, the bang-bang control switches values at t = tj\ j. 
Define 

The arguments in the preceding two paragraphs may be repeated until the initial 
time, t 0 is reached. Recall that tj~t 0 . Define v 0 = — )Vj and recall that previous 

A 

definitions require X x (r c )= these definitions imply satisfaction of (3.19). The 

proof of the “only if’ part is complete. 

For the “if’ pan of the theorem, assume that (3.66) satisfies {/}. Define 
t e and recall that and t } =t 0 . Define v, = -v 0 and v, = v ; . 

Given assumption (i), it is immediately obvious that all conditions in {///} except Eqs. 
(3.59), (3.62), (3.65), (3.61) with tel, and (3.60) with teN. Note that (3.61) and (3.60) 
each apply at a switching point and when u^u^. Furthermore, Eq. (3.14) specifies that 
the Hamiltonian is zero throughout the trajectory. Therefore, by Lemma ID. 1, Eqs (3.14), 
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(3.21), and assumptions (v) and (vii) there exists a different value v for each switching 
point such that Eqs. (3.60) and (3.61) hold; however. Lemma ID.1 does not guarantee that 
the value of v at one end of the *th «= 0 arc (/=*-l in (3.60)) equals the value of v at the 
other end (i=i in (3.61)). But. Lemma III.2 wi,hT(x(/)) = (»(,)) ,„ d assumption (vi) 

implies that X,(r)-^^~(x(/)J v solves (3.15) when u=0. Therefore, the value of v at 
one end of a u = 0 arc must equal the value of v at the other end of the u = 0 arc. 

Eq. (3.65) is implied by the switching function vanishing at t-tf. Finally, it is 
obvious that the boundary value problem cannot be solved if X t (t) = 0; therefore |v 0 | > 0, 

by assumption (v). That implies satisfaction of Eq (3.59). B 

ni.3.2. MPM Example Solutions 

The following examples satisfy all the conditions implied by the Euler-Lagrange 

equations and the Pontryagin Maximum Principle. All quantities have been 
nondimensionalized. 

The first example solution is a 5-bum transfer reproducing a solution presented in 
a paper by Redding. Both the initial orbit and the final orbit are circular. However, there 
is an inclination of 28.5° between them. In this presentation of the solution, the initial 
orbit is equatorial and the final orbit is inclined 28.5°. The initial orbit radius is 1, the 
final orbit radius is 6.4. The initial nondimensional acceleration is 0.0517 and the 
nondimensional characteristic velocity is 0.567. Both the transfer computed by Redding 
and this solution calculated with the Modified Patched Method have final transfer orbits 
with e =0.723 and an inclination 26.5° away from that of the final orbit. Perigee bum 
durations for both range from 1.26 to 1.13. Both have a total transfer time of 60. Finally, 
it is worth noting that the solution presented here was computed without knowing the 
particulars of Redding’s solution. 
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Orbit Number 

Components of the Angular Momentum Vector for each Transfer Orbit 
vs. Orbit Number of a 5-Bum Transfer with Plane Changes 



Orbit Number 

Figure 3.6 Components of the Eccentricity Vector Vs Orbit Number for each 
Transfer Orbit of a 5-Bum Transfer with Plane Changes 
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Figure 3.7 


Transfer Time Vs Burn Number 



Burn Number 


Numbcr for cach Burn of a 5 ' Bum Transfer 


The second example is a 19-burn transfer. The initial nondimensional 
acceleration produced by the rocket motor (T/m 0 ) is 0.09698 and the initial 
nondimensional characteristic velocity { go J sp ) is 0.3929. The initial orbit is circular with 
a radius of 1, the final orbit has eccentricity of 0.73315 and a semimajor axis of 9.26. 
The total bum time for this trajectory is 26.84. Figures 3.8 - 3.9 show data in similar 
form for this transfer as Figures 3.5-3.7 for the previous transfer. 


This 19-burn trajectory was extended to a 27-burn trajectory. This process 
involved the determination of transfers with 20. 22, 23, 24 bums, etc. It was found that 
adding bums one at a time was usually successful, two at a time slightly less successful, 
and so on. It was also interesting to see the decreasing improvement of the transfer’s 
performance as plotted in Figure 3.10. 
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Orbit Number 

Figure 3.8 Orbital Elements for each Transfer Orbit Vs Orbit Number of a 19-Burn 
Transfer 


Transfer Time Vs Burn Number 



* * ■ « » » * i • • • I i i l T ' I ' I ’ I 
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Burn Number 

Figure 3.9 Transfer Time vs Orbit Number for each Bum of a 19-Bum Transfer 
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Number of Burns 


Figure 3.10 Performance Index, Final Mass, of the Extremal Trajectory vs Number 
of Bums Executed during the Trajectory 

The third example is the aforementioned 27-bum trajectory. All parameters are 
identical between this transfer and the previous except the number of bums. The total 

bum time for this trajectory is 26.64. This is only a 0.7% decrease in transfer time for 
42% more bums. 



Orbit Number 

Figure 3.11 Orbital Elements for each Transfer Orbit Vs Orbit Number of a 27 
Bum-Transfer 
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Transfer Time Vs Burn Number 



Burn Number 

Figure 3.12 Transfer Time vs Orbit Number for each Bum of a 27-Bum Transfer 

The fourth transfer is identical to the third except that the final orbit has an 
inclination of 63.4° This inclination angle was chosen because it is large and represents 
the inclination of the useful Molniya class of orbits. To obtain the solution, the planar 
transfer was used as the initial guess and the Modified Patched Method obtained the 
solution in 6 iterations. The following figures represent the transfer. 

Each of these transfers show similar trends. An almost linear variation in the 
largest components of the angular momentum and eccentricity vectors and for the transfer 
time when plotted against the orbit or bum number. However, this trend is broken for the 
last bum. In each transfer, the last bum is an apogee bum and all previous bums are 
perigee burns. Each perigee burn steadily changes the angular momentum and 
eccentricity. The apogee bum then makes a last large change that brings the spacecraft to 
the final orbit. This last bum is also considerably longer than the bum before it. In the 5- 
bum case, Fig. (3.7) shows that the last bum is much longer than the first bum. In the 19- 
bum case, Fig. (3.9) shows the last burn almost just as long as the previous bum; in the 
27-bum case, Fig. (3.15) indicates that it is considerably longer. 
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Orbit Number 

Components of the Angular Momentum Vector for each Transfer Orbit 
vs Orbit Number of a 27-Bum Transfer with a 63.4° Plane Change 



Orbit Number 

Figure 3.14 Components of the Eccentricity Vector vs Orbit Number for each 
Transfer Orbit of a 27-Bum Transfer with a 63.4° Plane Change 
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Figure 3.15 Transfer Time vs Orbit Number for each Bum of a 27-Bum Transfer 
with a 63.4° Plane Change 

One feature that seems common to the large number of bums case and the small 
number of bums case is the use of the distant bum for inclination changes. Referring 
back to the nonplanar 3-bum transfer shown in Figs. 2.2-2.3, it is clear that the first bum 
is making most of the inclination change. Also, it is clear from the 27-bum transfer 
represented in Fig. (3.13) that the h y component of the angular momentum vector, which 
indicates the inclination, has very little variation until the final bum takes its value from 
almost zero to almost -2. This same trend can be seen for the 5-bum transfer represented 
by Fig. 3.5; where the h x component indicates inclination for this transfer. 

III.4. Inclusion of Pert urbation Terms 

Neither the Patched Method nor MPM are equipped to produce exact solutions to 
fuel-optimal orbit transfer problems in the presence of orbit perturbations. Note that 
including orbit perturbations will cause assumption (i) from {P} to be violated. 

The tradeoff between making the ideal gravity assumption and obtaining solutions 
w-ith much larger numbers of bums was deemed acceptable. It is hoped that the 
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techniques used in this tradeoff will find application in future research into the orbit 
transfer problem including perturbations. 

However, BOUNDSCO was able to obtain a solution including orbit perturbations 
for the 5-bum transfer presented above in Figure 3.5. Perturbations are considered for 
this trajectory as opposed to the others, because BOUNDSCO iterations did not converge 
for the others, even after several trials including initial guesses that were slightly 
perturbed from the exact solution. 

Figures 3.16-3.18 shows the changes in orbital elements and transfer time induced 
by the inclusion of atmospheric drag and oblateness effects. It is clear that the extremal 
trajectory includes a lengthened second bum which raises the energy of the second 
transfer orbit, thereby raising its altitude and decreasing the effect of drag. It is not so 
clear what decides that the longer bum will be the second and not the first. The nodal 
regression seems to manifest itself as a decreasing H y component; it is interesting to note 
that, like inclination changes, the extremal transfer doesn’t make the correction until the 
last bum. Turning attention to the bum lengths, note that the amount by which the first 
burn is shortened almost exactly counters the amount by which the next burn is 
lengthened. A similar trend shows itself for the third and fourth bums. The last bum is 
only slightly shorter, but not enough to indicate whether the total bum time is longer or 
shorter. In fact the final mass of the ideal gravity transfer was 3.762; for the transfer with 
perturbations it was 3.760. This is a performance loss of only 0.07%, a surprising result 
considering that the individual bum times change by as much as 1.6%. 


81 






D 

CD 

O 

0.01 

CD 

CJ 

V) 

CD 

0.008 

5‘ 

O 

0.006 

g 

sj 

0.004 

Io 

^3 

CD 

3 

0.002 

tn 

8 

3 

0 

5 

CD 

CD 

-0.002 

o 

CD 

in 

a 


0 1 2 3 4 5 g 

Orbit Number 

Decrease in the Components of the Angular Momentum Vector for the 
5-Bum Transfer of Figure 3.5 Considering Orbit Perturbations 





0.012 


0.008 


0.006 mm 
.* o 
m3 

0.004 -tta 


0.002 


- 0.002 


2 3 

Orbit Number 


Decrease in the Components of the Eccentricity Vector for the 5-Bum 
Transfer of Figure 3.5 Considering Orbit Perturbations 


82 




Figure 3.18 Decrease in the Bum Times for the 5-Bum Transfer of Figure 3.5 
Considering Orbit Perturbations 

HI.5 Conclusions 

In this section, two new methods for computing multiple-bum orbit transfers are 
presented. These methods, the Patched Method and the Modified Patched Method, have 
been developed specifically to fill an apparent gap in computational ability for fuel- 
optimal transfers with large numbers of bums. For this type of problem, both methods 
have out-performed BOUNDSCO and MBCM from the previous section. 

The conditions upon which each of these methods are based on have been proven 
equivalent to necessary conditions. However, for both methods it is required that 
Ponoyagin’s Maximum Principle be checked after iterations have stopped. 

The Patched Method, though slow, was very robust in obtaining solutions. 
Because of its use of a direct method, it was usually able to obtain the one-bum solutions 
between each pair of orbits. Also, the optimization of the transfer orbits usually 
proceeded well in the sense that each iteration would produce a better choice of orbital 
elements. However, the overall method tended to be quite slow because the cumulative 
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time required to compute the one-burn transfers in succession was quite long and 
increased with number of bums. 

MPM computed solutions beyond the capability of any of the other methods 
investigated in this report. MPM was much quicker and slightly less robust, as would be 
expected of a method more akin to multiple-point shooting. Therefore, it is suggested 

that the Patched Method be used with a very low tolerance to obtain initial guesses for 
MPM. 


Neither the Patched Method nor MPM is designed to handle orbit perturbations. 
However, the marked improvement in performance found with these configurations 
should be motivation enough for a future research effort to produce similar configurations 
that can handle orbit perturbations efficiently. 

Also in this section, a new formulation for the solution of the Lagrange 

multipliers is presented. This formulation is valid over coast arcs where the Hamiltonian 
vanishes. 


84 



Section iv 

Guidance for Optimal Orbit Transfers 

IV.l. Introdnrfinn 

The guidance scheme examined here is an implicit one which implements 
neighboring optimal feedback guidance. An implicit guidance system was chosen due to 
the fact that that type of guidance system often handles disturbances well 42 . Neighboring 
optimal feedback guidance was chosen because it has the advantage of being a feedback 
system, as opposed to open-loop guidance and it can be implemented very easily as with 
a gain-scheduling scheme. There also appears to be a lack of studies in the literature 
which examine this type of guidance scheme for this problem. 

In this formulation, the initial orbit exit point is assumed to be perturbed from the 
nominal point but the other boundary condition, specifying the final orbit, is assumed 
unchanged. The goal is to use the controller to bring the trajectory to the final orbit at 
some point with minimal fuel. 

In order for this guidance scheme to be implementable, the neighboring trajectory 
must exist; the sufficient conditions for a local extremal must be satisfied. The 
satisfaction of these conditions for the nominal solution will be shown. Following that, 

the guidance scheme will be investigated, including the use of a time-to-go indexing 
scheme. 


42 Naidu, D. Subbaram. Aeroassisted Orbital Transfer: Guidance 
New York: Springer-Verlag, 1994. 


and Control 


Strategies. 
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IY.2. Literature Rfyjfty 

Many researchers have used the first variation to compute extremal solutions to 
the fuel -optimal orbit transfer problem. However, few, if any, have made use of the 
conditions related to the second variation of the cost functional in computation. These 

provide sufficient conditions which, when met, declare an extremal solution as a locally 
weak optima] solution. 

Once the second variation of the cost functional is verified so that it is known 

whether the sufficient conditions are met, the information obtained can then be used to 

implement a guidance scheme. Guidance schemes can typically be divided into two 

categories: implicit and explicit. Implicit guidance systems are characterized by the fact 

that the vehicle’s motion must be precomputed on the ground and then compared to the 

actual motion. The equations which need to be solved are based upon the difference 

between these measured and precomputed values. The solutions to these equations are 

used in the vehicle’s steering and velocity control. Explicit guidance systems are 

generalized by the fact that the vehicle’s equations of motion are modeled and solved for 

by on-board computers during its motion. The solutions for the equations are solved 

continuously and are used to determine the difference between the vehicle’s current 

motion and its destination. Commands are then generated to alleviate the anticipated 
error. 


Guidance schemes have been presented in various papers. 43 A guidance scheme 
which is implemented using a linear tangent law is presented by Sinha, Shrivastave, Bhat, 


43 Chuang C.-H., Goodson, TX>. Ledsingcr, L.A., “The Second Variation and 

^A?K»0^/ b /a k r Gu i danCe for ^ u - ltiplc Bum Orbit Transfers,” 

BatoOT. USA Cmfere " Ce ° n Guidmce ’ Navt * mim ' <”•“ C °™°'- 
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and Prabhu. 44 In a paper by Lu 45 * a nonlinear guidance law is developed using two 
different strategies. One strategy uses optimal control theory to generate a new optimal 
trajectory onboard from the start, while the other uses flight-path-restoring-guidance to 
bring the trajectory back to the nominal. A guidance scheme that is developed using 
inverse methods for unthrusted, lift-modulated vehicles along an optimal space curve is 
presented by Hough. 4 * Linearized guidance laws applicable to many different types of 
space missions are presented by Tempelman 47 These guidance laws are based on fixed 
and free final time arrivals. Naidu 42 presents a neighboring optimal guidance scheme 
applicable to aeroassisted orbital transfers. 

m. Preliminary Considerations 

Earlier, the optimal orbit transfer problem was given as a maximization problem. 
To conform to the convention used for the second variation 36 it is transformed to a 
minimization problem. For the minimization problem, the performance index can be 
made negative and considered a cost functional 


J = 



(4.1) 


As the necessary conditions are First-order conditions, they remain unchanged. However, 
Lawden s pointer vector theory is second-order and requires that the control be such that 


«r = 


-X 

\K 


(4.2) 


^Sinha, S. K., S. K. Shrivastava, M. S. Bhat, and K. S. Prabhu. “Optimal Explicit 

ro U 0 o CC £c 1 ??f' Dimensional Launch Trajectory,” Acta Astronautica. Vol. 9, 
iysy, pp. H5-125. 

45Lu ’ P a/ A - Ge . neral Nonlinear Guidance Law,” Proceedings of the the AIAA Guidance, 
Navigation, and Control Conference , Scottsdale, Arizona, 1994. 

46 Hough, M. E., Explicit Guidance Along an Optimal Space Curve,” Journal of 
Guidance, Control, and Dynamics. Vol. 12, 1989, pp. 495-504. 

47 Tempelman, W., “Linear Guidance Laws for Space Missions,” Journal of Guidance 
Control, and Dynamics. Vol. 9, 1986, pp. 495-502. 
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Furthermore, Pontryagin's Minimum Principle requires that an extremal solution satisfy 


where 


^<o, r = r_ 

H s > 0, 7 = 0 


(4.3) 


"i — 




m 



(4.4) 


If an extremal solution to the maximization problem is given as state time history 
x(t). La grange-multiplier time history A.(t), and Lagrange multipliers v, (associated with 

boundary conditions) then an extremal solution for the minimization problem with the 
cost function in Eq. (4.1) can be constructed as x(r), (-l)*L(r), and (-l)*v. 

Additionally, it makes more sense in the planar guidance problem to consider the 
control as an angle 6 , rather than individual components of a unit vector. This simplifies 
analysis because the control is now a scalar. Equation (4.2) now gives 


tan(6) = -j± (4.5) 

A practical approach to guidance is suggested by previous results in this report. If 
a multiple-bum transfer can be thought of as consisting of multiple optimal one-bum 
transfers, then it should be reasonable to examine a guidance scheme that attempts to 
match each of the intermediate transfer orbits of the multiple-burn transfer. In other 
words, use neighboring optimal feedback guidance for one bum at a time. 
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This is not suggested to be an optimal guidance scheme. By focusing on each 
bum with neighboring optimal feedback, but not considering the trajectory as a whole, 
this guidance scheme becomes a sub-optimal guidance scheme. 

Each bum can be considered an extremal solution. These extremal solutions are 
considered to have a fixed initial point and free transfer time but the final point is only 
constrained in that it must lie on the final orbit. Recall, however, that in computing the 
multiple-bum transfer the initial point was not fixed; this condition is imposed for 
practical considerations. If the spacecraft is delivered to the correct orbit, and coasting to 
the nominal bum-on point has zero cost, then there is no reason to attempt to compute a 
new bum-on point. This reasoning holds for the beginning of each bum. 

IY.4. The Second Variation for Ono -Burn Prohlprm 
Considering the second variation of the augmented cost functional, J, a new 
optimal control problem can be stated. 36 In this new problem, the state is <5x, the control 
Su, and the Lagrange-multipliers are Sk and dv. The new problem is linear and can be 

solved using a sweepback method. For the problem considered here, x=[r T v* m jT and 
u=6. 


When the final time is free for optimization, the transversality condition must be 
satisfied by the nominal solution. The notation for this condition is 


where 


n(*,v,i)| =f— 1 =(— *1 =o 


(4.6a) 


G(x, v) = <}>(x) + v T \y(x) 


(4.6b) 
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In general, neighboring optimal feedback guidance allows consideration of 
changes in boundary conditions. No such changes are considered, assuming that the 
destination orbit is fixed. Formulation will be made below for the free final time case. 


The change in state and costate can be estimated with a linear time-varying 
dynamic system. This dynamic system is given below, where it is understood that matrix 
functions arc evaluated with the nominal trajectory. 


^■<5x = A(r)<5x-B(r)<5X 

(4.7) 

~^ = -C(t)5x-A r (t)5X 

(4.8) 

> 

II 

M 

1 

h 

(4.9) 

BW'fA'fl 

(4.10) 

C(0 = H a - 

(4.11) 


Evaluating Eqs. (4.7)-(4.U) the recurring terms in the differentia] equations are: 


0 

0 

-(?)-¥ 

3/ay 
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-|T 

0 

0 

0 

3fjxy' 
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1 

1 

|u> 
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0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

cos(6) 

— — sin(6) 

0 

- 


m 

m 

- 

r _ 
e “ 

[» 

0 — — sin{6 ) 
m 

T 1 T 

— cos(8) 0 
m J 



(4.12) 


(4.13) 
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(4.14) 




Q o 
0 2 


Q = 2if.f{( 3A - Jr + A - >V - 5(>-:r> ! } {(A„y + A,*)r ! - 5(Ajr)a>] 

[ {(A_v + A.x)r 3 - 5(XTr)^} {(3A.y + A.x)r , -5(tf r )jry}J 


w„=-|x, 

m 

=0 


(4.15) 

(4.16) 

(4.17) 


note that r = [ar y]\ v = [„ vf, and X. = [A. A.] 1 are taken as the nominal 

trajectory. Using the sweepback method for nonlinear terminal constraints the form for 
57, and 5y can be written as 


5X(r) = P (t)5\{t) + S(t)dv 

(4.18) 

Sy = S T (r)<5x(r) + V(i)dv 

(4.19) 

which allows the solution for dv to be written as 


dv = V"’(r,)[5t|/ - S T (r,)$x(r,)] 

(4.20) 

As mentioned above. S\y=0 will be considered here. The matrices F(t), 
are computed using the following relations: 

S(/), and V(r), 

P(r)-P(r)- n,( ' )n,T W 

a{t) 

(4.21) 

vj| 

*-» 

II 

on 

i 

3 

S C, 

(4.22) 

V(t) - V(t) "«" T W 
a(t) 

(4.23) 
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Now the matrices P(r), S(t), V(r), m(r), n(r), and the scalar function a(r) are computed 
from a dynamic system. The boundary condition equations for this system are given by; 


p ('/) = [<'= + (v T V,) I ],., < 
s ('/)=[vl],., ; 
v(,,)=o 

where in the development for the orbital transfer these are; 

a b d e O' 

b c f g 0 

d f h i 0 

eg i j 0 

0 0 0 0 0 




d = -v 3 v 

e - Vj - v 2 u + 2 v 3 v 
f = -v l - v 2 v + 2 v 3 u 


g = ~v 2 u 

h = 2v 3 y 

/ = -v 3 ;c- v 2 y 
j = 2 v 2 x 

and expression for Eq. (4.25) was previously given as Eq. (2.1 1). 



(4.24) 

(4.25) 

(4.26) 


(4.27) 


(4.28a) 

(4.28b) 

(4.28c) 

(4.28d) 

(4.28e) 

(4.280 

(4.28g) 

(4.28h) 

(4.28i) 

(4.28J) 
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Following from the assumptions expressed as Eqs. (4.18)-(4.19), the following 
nonlinear equations for P, S, and V must be integrated backwards. The results will be 
used to check the sufficient conditions governing a minimizing solution. 


-PA- A t P + PBP-C 

(4.29) 

s = -(a t -pb)s 

(4.30) 

v=s t bs 

(4.31) 

m = -(A T -PB)m 

(4.32) 

n = S T Bm 

(4.33) 

a = m T Bm 

(4.34) 


with the following boundary conditions applying 


"M-(f 

) 

(4.35) 

■«■(?) 

~v 

'*'/ 

(4.36) 


'"i 

(4.37) 

The sufficient conditions for a minimizing solution can now be stated as follows: 
convexity condition: H ee {t) >0 for t,<>t£t f 

(4.38) 

. V -1 (r) exists for t 0 <t<t, 

normality condition: 1 

(4.39a) 

a‘ (t) exists for t 0 < t < t J 


(4.39b) 

conjugate point condition: P(r)- S(r)V'(/)S T (0 finite for 

(4.40) 
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The convexity condition is satisfied for any transfer satisfying Equation (4.5). This can 
be seen by noting that Eq. (4.16) is positive definite, irrespective of the time history for 
the Lagrange multipliers. 

The eigenvalues of V are plotted in Figure 4.1. Figures 4. 2-4.4 plot the elements 
of the conjugate point condition matrix. Figure 4.5 is a plot of a(t). Figure 4.1 shows 
that V is positive definite in the required interval. Figure 4.5 shows that a(t) is negative 
definite in the required interval. Since the normality condition requires that the inverse of 
V and a(t) exists in the interval, this solution is normal. Figures 4. 2-4.4 show that the 
conjugate point condition is satisfied. The elements are bounded in the required interval 
and grow asymptotically at the final time; the curves in the figures have been truncated to 
show their variations prior to this asymptotic growth. Therefore, this solution satisfies 
the sufficient conditions for minimizing the cost functional with free transfer time. 

It seems appropriate to first attempt the guidance scheme for a relatively 
uncomplicated transfer. Such a transfer was presented in Fig. 2.1 and discussed in 
subsection [11.2.4], The transfer is planar; no plane changes occur. The guidance scheme 
considered here will be simulated for this trajectory. 


IV.4.1. Neighboring Optimal Feedback Guidance 

Conveniently, construction of a neighboring optimal feedback guidance law uses 
the same information as that required to check the second variation of the cost functional. 
As a result, much of the derivation required of guidance law has been stated already. The 
remaining discussion will describe how to form the feedback control law and adjust the 
characteristics of the bang-bang control in a feedback law. 

The control, 66, for the fixed final time problem can be found using 
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Figure 4.1 



time (normalized) 


Plot of Eigenvalues of V(r) for Two-Bum Extremal, Last Bum 



time (normalized) 

Figure 4.2 Plot of Elements of Conjugate Point Condition Matrix for Two-Bum 
Extremal, Last Bum 
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o 0.2 0.4 0.6 0.8 1 

time (normalized) 

Figure 4.3 Plot of Elements of Conjugate Point Condition Matrix for Two Bum 
Extremal, Last Bum 



o 0.2 0.4 0.6 0.8 1 

time (normalized) 

Figure 4.4 Plot of Elements of Conjugate Point Condition Matrix for Two Bum 
Extremal, Last Bum 
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Plot of a(t) for Two Bum Extremal, Last Bum 

This continuous feedback law has been constructed by estimating dv at each instant of 
time instead of solving for dv at the initial time and then using this value for all time 

The feedback law depends on P, S, and V as functions of time. A particular 
advantage of neighbonng optimal feedback is that the linearized TPBVP only has to be 
solved once. Afterwards, sampled values of the feedback gains may be stored. The 
feedback gains may then be computed for any time by interpolation between stored 

values. Use of this control should keep the spacecraft on a neighboring optimal solution 
and deliver it to the required orbit. 

The block diagram for the feedback controller needed for neighboring optimal feedback 
guidance is shown in Figure 4.6. 
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Figure 4.6 Diagram of Neighboring Optimal Feedback Controller Implementation 


where Ai(t) is the feedback gain from Eq. (4.41), computing 50. 

IV.4.2. Simulation of the Guidance Algorithm 

Justification for a feedback algorithm lies in Fig. 4.7 and Fig. 4.8. It can be noted 
that there is error in the variation of the states from the neighboring optimal trajectory 
when guidance is not used, Fig. 4.7, i.e., when the control correction is not used. 
However, Fig. 4.8 shows that a feedback law is needed because when implementing it, 
the errors in the variation of the states becomes much less, comparatively, than that using 
no guidance whatsoever. The neighboring optimal trajectory referenced in Figs. 4.7-4.8 
was computed with BOUNDSCO. 

IV .4.3. Time-To-Go Implementation 

Since this problem is a free final-time problem, the possibility exists that the final- 
time will increase and the guidance algorithm will “run out of gains”; this is a familiar 
issue for neighboring optimal feedback guidance. The approach used in this study is 

based on discretizing the gains by N time nodes where t N is earlier than the 

nominal if The gains at the nominal ^ will be infinite and impractical to store. Both the 
gains for calculating dt } , via Eq. (4.42), and for 60, via Eq. (4.41), are then calculated at 
any time by linear interpolation between stored values. - 


98 





Plot of State Variation from the Nominal Trajectory vs. Time for 
Neighboring Optimal Trajectory ((),*„) and a Trajectory Without 
Guidance (no subscript) 


Figure 4.8 



0 0.2 0.4 0.6 0.8 

Time (normalized) 


Plot of State Variation from the Nominal Trajectory vs. Time for 
(no subscript) Trajectory (() ncn ) and a Trajectory With Guidance 
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To consider time-to-go, the guidance must make active use of the dtj estimation. 
Since both the nominal and the actual trajectories stan at t lt dt fl can be initially calculated 
using the gains at that time. The length of the first guidance interval is then found by 
relating it to the estimated time-to-go. 


t.+di , . 

Ar„ = H h~h) 


(4.43) 


Then, at the end of the /-7th guidance interval, the gains at /, are used to calculate di A . 
Using this information, the length of the ith guidance interval can be computed as 


i-1 


'/X-2X 


Ar.. = • 


i* i 






(4.44) 


This continues until At, is computed as zero or a negative number or until N. When 
i=N, the Nth gain is used for the entire interval A t$. When this interval ends, the 
guidance scheme is finished. 


The plots below compare guidance performance with and without this time-to-go 
formulation. The curves represent the time history of the boundary' condition error, i.e. 
Eqs. (1.12) minus the desired orbital elements, evaluated continuously. Figure 4.9 makes 
continuous use of the gains but indexes these gains at the current actual time without 
calculating dt f . For the perturbation simulated, the transfer time needs to increase and this 
first scheme must terminate prematurely. Figure 4.10 makes use the discretized gains and 
time-to-go formulation. This simulation also incorporates a practical saturation limit on 
the size of the gains. The improvement due to the time-to-go formulation is obvious 
when comparing these plots. Therefore, this is both a practical and superior 
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implementation of the continuous bum guidance considering the boundary condition 
error. 



time (nondimensional) 

Figure 4.9 Plot of Boundary Condition Error for Continuous Guidance 



time (nondimensional) 

Figure 4.10 Plot of Boundary Condition Error for Discrete Guidance with Time-to- 
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1Y.5. Multiple Burn Guidance 

The guidance for multiple bums can also be discretized. For the two bum case, 
discretized guidance using time-to-go is used for the first bum. The guidance algorithm 
will place the spacecraft on the intermediate transfer orbit via the neighboring optimal 
trajectory. Since the cost on this coast arc is zero, the spacecraft can coast on this arc 
until it reaches the point at which the next bum is to start. Once the spacecraft reaches 
this point, discretized guidance using time-to-go can be used again for the second bum. 
The boundary conditions for the second bum should than be satisfied by the neighboring 
path. For multiple bums, this guidance scheme is extended in a straightforward manner. 

The guidance scheme detailed above was used to recover the two bum transfer of 
Fig. 2.1 in the presence of an initial perturbation. Fig. 4.1 1 shows the boundary condition 
errors for the first bum given an initial perturbation of 10'^ in non-dimensionalized units. 
The boundary conditions are satisfied rather well for this bum. The resulting boundary- 
condition errors for the second bum are shown in Figure 4.12. The boundary conditions 
are satisfied very well for this bum. 

Figures 4.13 & 4.14 show the boundary condition errors during the second bum 
for a perturbation of the same magnitude as above in only the x position and the u 
velocity, respectively. Note that the error in the boundary conditions is slightly greater in 
Figure 4.14. This suggests that the trajectory is more sensitive to disturbances in the u 
velocity than in the x position. 
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0 1 2 3 4 5 6 7 8 

time (nondimensional) 


Figure 4.11 Plot of Boundary Condition Error for Discrete Guidance During the 
First Bum 



Figure 4.12 Plot of Boundary Condition Error for Discrete Guidance During the 
Second Bum (Continuation of Fig. 4.11) 
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time (nondimensional) 

Figure 4.13 Plot of Boundary Condition Error for Discrete Guidance During the 
Second Bum for error in x at the initial time 



time (nondimensional) 

Figure 4.14 Plot of Boundary Condition Error for Discrete Guidance During the 
Second Bum for error in u 


The resulting orbit transfer trajectory is shown in Figure 4.15. This plot corresponds to 
the boundary condition errors as shown in Figures 4.1 1 and 4.12. 
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X (km) 


Figure 4.15 


Plot of the Two Bum Orbit Transfer (from Fig. 2.1) with Initial 
Perturbation 


Extremal one bum trajectories have been shown to be weak locally optimal 
solutions using sufficient conditions. This docs not prove that the multiple-bum transfer 
from which they were taken is itself a weak locally optimal solution, but it does allow the 
use of a new subopdmal guidance scheme. 

This scheme was shown to reduce the terminal errors for small perturbations of 
the initial state. To increase the size of allowable perturbations, a time-to-go indexing 
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scheme was simulated. This time-to-go indexing did improve the performance of the 
guidance scheme. 

The suboptimal multiple-bum guidance with time-to-go indexing was simulated 
ior a planar transfer. The performance of this guidance scheme did not match 
expectations. The implication is that the region in which a linear control correction is a 
valid assumption was quite small. Actually, this is not a surprising conclusion since 
obtaining the nominal solutions is usually quite a challenge for iterative algorithms that 
attempt linear corrections for each iteration. If indeed this implication is correct, then a 
more sophisticated approach for neighboring feedback control is required. 
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Section v 


Conclusions and Recommendations for 
further Study 

y.l. T ransfers With Small Numbers of Burns 
It has been found that methods already present in the literature are capable of 
computing fuel-optimal orbit transfers with small numbers of bums. The methods 
investigated here were multiple-point shooting and modified shooting. However, a 
common way to attempt to increase the performance of a transfer is to increase the 

number of bums executed and, unfortunately, these methods are not very robust in that 
sense. 


A new method has been introduced that is very useful for adding bums to fuel- 
optimal orbit transfers. The method is used in conjunction with homotopy and an 
iterative technique for computing transfers; the iterative technique must incorporate 
knowledge of the Lagrange multipliers. The method does require that the initial point, 
the final point, and the transfer time be free for optimization. It also assumes that the 
transfer is performed under the influence of ideal gravity. This assumption is required to 
obtain the switching function property that the method relies on. 

It is recommended that this method be further developed such that orbit 
perturbations are taken into account. Since the switching function property in question 
no longer applies for this case, the task is challenging. Obviously, a fairly different 
approach must be taken. It is likely that requiring trajectories to begin and end with coast 
arcs will be necessary, since cost arcs will no longer be orbits. Perhaps then some 
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conditions may be identified under which the coast arcs could be extended to find optimal 
locations for the bums to be added. 

Transfers with Large Numbers of Burns 

The results of this research point to the Modified Patched Method as a practical 
way to compute fuel-optimal transfers with large numbers of bums. It does not appear 
that such a method existed previously in the literature, making MPM and theoretical 
results behind it the central contributions of this report. 

An interesting spin-off of the theoretical development is a new formulation for the 
integration of the Lagrange multipliers over a time-optimal coast arc for the nonplanar 
case assuming ideal gravity. The formulation results from satisfaction of Lemma III. 2. 
This particular formulation proved Quite useful for MPM and may prove useful in future 
algorithms and future theoretical developments. 

MPM does not allow for orbit perturbations. This restriction was a small price to 
pay for performance previously unobtained, viz. the ability to compute transfers with 
upwards of 27-bums and large inclination changes. Now that this performance has been 
obtained for the ideal gravity case, it is suggested that a future research effort should be 
able to produce a method with similar performance, or better, while taking orbit 
perturbations into account. 

If an attempt is made to adapt MPM for orbit perturbations without recovering 
any properties lost, then MPM will degenerate into multiple-point shooting. This study 
has already concluded that multiple-point shooting does not perform well for large 
numbers of bums; therefore, some recover)' of the properties from Theorem ni.l and/or 
Theorem 111.2 must be made. Since the concept central to both the Patched Method and 
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N1PM is the relationship of the optimal orbit transfer problem with the problem expressed 

by (3.1), it seems reasonable to expect some form of (3.1) to be recovered in the presence 
of orbit perturbations. 

¥.2. Multiple-Burn fimdanrp 

A suboptimal multiple-bum guidance scheme was developed through this 
research and its performance investigated. The scheme may be described as "bum-by- 
bum" neighboring optimal feedback guidance with a time-to-go indexing scheme for each 
bum. The performance of this guidance scheme did not match expectations. 

Since guidance has much practical importance, it is suggested that future research 
attempt to develop an improved guidance scheme. It is likely that this would involve 
techniques to improve neighboring optimal feedback or replacing this with some other 
one-bum guidance scheme. On the other hand, a future research effort might attempt to 
find an optimal guidance algorithm for the multiple-bum transfer as a whole. Since there 
is a strong relationship between the sufficient conditions for optimality and the 
computation of neighboring optimal feedback gains for the one-bum problem, a similar 
relationship might be expected for the multiple-bum problem. If an optimal multiple- 
bum guidance scheme is developed, it will likely lead to the development of sufficient 
conditions for the optimality of multiple-bum transfers. 
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ORB PACK Users Manual 


A Package of FORTRAN Programs to Construct 
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Applied Control Laboratory 
Georgia Institute of Technology 
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L Introduction 

ORBPACK is a collection of FORTRAN 77 programs for computing 
optimal orbit transfers. For the most part, these are all indirect 
methods; they are concerned with solving the Two Point Boundary 
Value Problem provided by optima] control theory. 

None of these routines guarantee a globally optimum solution, only 
extremal solutions are claimed by convergence of iterations With 
the exception of MBCM, solutions obtained with these methods must 
have their switching law checked. One must be sure that, in the 
computed solution, the thrust is on when the switching function is 
positive and the thrust is off when the switching function is 
negative. Furthermore, these methods assume that no intermediate 
thrust arcs will be found in the solution. 

The charts below summarizes the programs in ORBPACK 


Solvers 



Method 

Libraries 

Suggested Use 

BND3D 

Multiple Shooting 
(BNDSCO) 

BNDSCO 

medium/low thrust; 
few burns 

MB CM 3D 

Shooting w/ Minimizing 
Boundary Condition Method 

VF02AD 

medium/low thrust; 
few burns 

PAT2D 

Patched Method 

BNDSCO; IMSL 

medium/low thrust 

MPMM2D, 

MPMM3D 

Modified Patched Method 

IMSL; ODEPACK 

medium/low thrust; 
short burns 


Accessories 



"Use 

Libraries 

rtirm— 

random shooting for one-burn guesses 

IMSL; ODEPACK 

p? 1 1 1 1 1 

convert MPMM2D files to MPMM3D files 

N/A ~ ' 


[convert MPMM3D files to BND3D files 

ODEPACK 1 

| BXU2MBCM 

convert BND3D files to MBCM files ' 

N/A 


All codes as supplied in ORBPACK solve multiple bum orbit 
transfers with free final time and free initial and final points 
BND3D is already configured so to switch between free and fixed 
final time problems. MBCM3D can easily be reconfigured for such. 
PAT2D, MPMM2D, and MPMM3D have fixed configurations. 

PAT2D, MPMM2D, and MPMM3D are also fixed to solve only 
problems where ideal gravity is assumed. BND3D and MBCM3D 
are configured to solve problems that include drag and oblateness 
effects. Finally, codes with the “2D” suffix are configured to solve 
planar transfers; the “3D” suffix indicates that the code is 
configured for nonplanar transfers. 
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H Orbit Transfer Problem Definition 


n.l. Parameters 


All the programs in ORBPACK require the following orbit transfer 
parameters to be determined 


For the gravitating body: 

• the gravitational constant for the central body (p ) 

For the rocket motor: 

• maximum thrust 

• specific impulse (/„„) 

sp 

For the terminal orbits, BND3D and MBCM3D require 

• semimajor axis 

• eccentricity 

• right ascension (degrees) 

• argument of perigee (degrees) 

• inclination (degrees) 


For the terminal orbits, MPMM2D, MPMM3D, and PAT2D require 

• angular momentum vector (X, Y, Z components) 

• eccentricity vector (X, Y components) 


tach program also requires a value for Earth’s acceleration at sea- 
level (g Q ) m appropriate units; this number is only used in 
conjunction with the specific impulse to compute the fuel 
consumption. 


BND3D and MBCM3D can account for oblateness and drag effects 
For oblateness: R e is the equatorial radius of the central body and J 2 
is a constant describing the mass distribution of the central body 
for Earth J 2 =1082.61xl0-6. For drag: fi is a constant from the ’ 
atmosphere model describing air density variation in the 
prescribed altitude region, p 0 is the atmosphere density at the 
altitude Tq, S is the cross-sectional area of the craft, and C n is the 
craft’s drag coefficient. 


The gravitational potential, including oblateness, is modeled 

where r is the magnitude of the position vector r. The drag force 
modeled as: 6 e 


as: 


is 




where v is the magnitude of the velocity v. 
density variation indicates an isothermal 


Note that this form for the 
region of the atmosphere. 
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D*2« Scaling It is very useful for numerical methods to work with numbers that 

are at or near the same order. This can be accomplished through 
nondimensionalizations. Such nondimensionalizations for the 
orbit transfer problem follow: 

f - r/r* 

m m m/m * 

t m tf 

and they require the following 


V- v/V/i/r* 

'*/ 

■ rjr ° 


(A,5C D ).p fl 5C D (r7 W ^) 

r-(r/«*)/(Ai/r a2 ) 


(iJ v ) ■ gJsp^r* hi 

K • 

Note that these nondimensionalizations result in dynamics with 
/**!• The choices of and are completely arbitrary A choice 
for might be one such that the initial nondimensionalized mass 
is 1 or 10. A choice for r* might be the radius of the planet or a 
number such that the initial semimajor axis, radius of perigee , or 
an “average” radius is 1. 
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m. Making Guesses for the Optimal Transfer 

There are many different ways that one could conceive of to make 

pvjesses The routines for making guesses, listed below, have been 
provided. 


The tutorials in Chapter VIII demonstrate how to make guesses with 
these methods. 


m.l. GSHOOT Random 
Guess (Single Bum Only) 


The subroutine GSHOOT will randomly make guesses for the one- 
urn orbit transfer problem in two dimensions. Input for GSHOOT 
is a text file. Its output consists of two text files which represent data 
tor direct and indirect methods. 


How to use GSHOOT 


GSHOOT requires a file, named “GINPUT,” for input 
“GINPUT” file follows: 


A typical 


MU * 1.0C 

GO *1.00 

ISP * 0.5673 

THRUST * 0.5166 

MO * 10. COCO 

AC - 1 DC 0 CC- 

EC S C.000 

wo * 0.000 

AC . 1.285 

ED * 0.219 

WD * o.ooc 

TKAX = G . 000 

NGS « 100 

NIX ■ 3 


where MU (n ) is the gravitational constant, GO (g 0 ) is the 
gravitational acceleration of the earth at sea level ,°ISP (7 ) is the 
motor’s specific impulse, and Thrust is the motor’s thrusHevel MO 
(m 0 ) is the initial mass for the transfer. The next parameters 
specify the terminal orbits: AO ( a 0 ) is the initial orbit’s semimajor 
axis, EO (e 0 ) the initial orbit’s eccentricity, and WO (w ) is the 
initial orbit’s argument of perigee; AF (a f ), EF (e f ), and°(ay) are the 
corresponding parameters for the final orbit. TMAX is the 

rJSnnT tim6i if 4t j 8 SCt to zero ’ then TMAX is assigned by 

„^ 00 . T t0 the amount of time required for the mass to vanish 
NGS is how many guesses to make; half of these will be almost 
tangential thrusting with random initial true anomaly and the 
other half will have random initial direction and random initial 
true anomaly. For a detailed description of the file format, see 
Appendix A. 


W1 Create out P ut fi,es “DIRECT.DAT” and 
i IRECT.DAT” which can be used to construct a multiple burn 
guess in the PATCH2D file format. Both of these files have 
identical headers: 



How GSHOOT 
works 


T * « 

GO . « 

ISP « f 

AC > i 

nxc « • 
DC • « 
AT * « 
EXT * « 
FYF « « 


These output files contain the necessary information 

If this output file represents a guess for any but the last bum delete 
the last three of these lines (AF, EXF, EYF) when constructing the 
multiple-burn guess file However, if this guess is for the last burn 

pvni rVf 1 three lmeS and delete Unes six thf ough eight (AO, EXO 

lines (T GO, ^SpV 5 ^ ^ ^ ^ delet€ the first thr «’ 

GSHOOT makes a random guess by choosing the constant 
agrange multipliers (v) as a random vector with unity magnitude 
Since all the Lagrange multipliers may be scaled by an arbitrary 
constant, there is no loss of generality. The state vector is computed 
knowing the initial orbital elements and randomly choosing the 
initial true anomaly. Next, the vectors ? r and ^ are calculated for 
the initial time, using the following equation: 



(3 1) 
v 


The initial value for A m is found by specifying that the switching 
function is zero at the initial time: 


m v„) 


(3.2i 


That the switching function is zero at the initial time is known to be 

the fre f transfer time and free terminal points problem. 

, ,he ‘ n ! Ilal Et ? te and costale known, the initial value problem i> 

1”, l 1° I"”' , ““ il ei,her the desired semimajor 

axis (AD) is reached, the current radius becomes small the 

spacecraft enters a parabolic orbit, or the mass becomes small 

For guesses that are almost tangential, \ is chosen to be (+/-) v and 
k is chosen to be (+/-) (ji/r 3 )r. The positive sign usually produces 
orbit raising and the negative sign orbit lowering Note that this 
initial guess for the costates zeros the Hamiltonian when the 
sw itching function is zero. Therefore, the v,’s can be found by 
solving the least-squares problem of Eq ( 3 . 1 ). 

Wil1 tJ T as many guesses as the user requests. The guess 
that best meets the required boundary conditions will be output 
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UL2. PAT2D Sub-Optimal 
Transfer Guess (Multiple 
Burn Only) 


PAT2D creates sub-optimal trajectories in the sense that the choice 
of intermediate transfer orbits has been fixed and each burn is an 
optimal one-burn orbit transfer. PAT2D iterates upon the choice of 
intermediate transfer orbits until it finds a choice that gives a local 
maximum in final mass The PAT2D program is described in 
detail in Chapter V. 


Usinp PAT2D to 
Compute Guesses 


P AT2D requires two files for input. The first file, 
“PATCH2D.TOLS,” sets accuracy levels and limits the number of 
iterations (for more information on this file, see Chapter V). The 
second file, “PATCH2D. GUESS,” supplies the guess information 
for both the choice of intermediate transfer orbits and the trajectories 
of the burn arcs between them. This latter file must be in the PAT2D 
format (for more information, see Appendix A and Chapter V). 


The guess information from GSHOOT, or some other source, must 
be put into the PAT2D format. When run, the first thing that PAT2D 
will do is solve the one-burn problems defined by the intermediate 
transfer orbits. Often, the output from this step alone is a 

“PATCH 20*1 NmAL 0 ” 1 111688 ^ ° utput is contained in the file 


On the other hand, it is not uncommon for that output to be an 
insufficient guess. In this case, one approach is to allow PAT2D to 

“PATruon^pcT^ during the iteration, the user may take the file 
rAiCH2U.bhbT and use it as an initial solution guess. 

the user ma >' set a rather loose stopping criterion for 

-PATrHon 5“’! “^"7°" is " e * In this approach, the file 
-rAiCxi2D.SOL will be the solution guess. 
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IV. The Modified Patched Method (MPMM2D, MPMM3D) 

The subroutine MPM2D (MPM3D) is a realization of the Modified 
Patched Method in two (three) dimensions The file “MPMM2D.T 
(MPMM3D.fi) contains an implementation of MPM2D (MPM3D; 
using IMSL s NEQNF to solve the nonlinear equations, its 
FORTRAN program name is MPMM2D (MPMM3D). ’ 


IV. 1. Using MPMM2D to 
Cc^i^ute Solutions 


MPMM2D (MPMM3D) requires only one input file, which must 
follow the PAT2D (PAT3D) format (6ee Appendix A). This data file 
must be named “MPM2D. GUESS” (“MPM3D. GUESS") 


The code “MPM2D3D.F will convert an “MPM2D. GUESS” file into 
a “MPM3D. GUESS” file. In this code, no other input is required 
except “MPM2D. GUESS” 


Data File (Input) In “MPM2D. GUESS,” (“MPM3D. GUESS”) the tolerance setting 

(TOL) is the root-finding tolerance. The tolerance used in 
numerical integration is one-thousandth of this number. No 
information in the header is ignored. 

For MPMM2D (MPMM3D), the option SEL may only be chosen as 1 
or 2. These options indicate the data for the bum is given in the 
format for an indirect method. MPMM2D (MPMM3D) will treat 
both SEL=1 and SEL=2 identically. 


MPMM2D (MPMM3D) only uses specific items from the PAT2D file 
format. The lines below are representative of the data for one burn 
m the PAT2D format. The underlined symbols indicate which 
number items are important to MPM2D calculations. 


MPMM2D 


a = 

i 






ex = 

1 






ey 

i 






NODE = 

3 






SEL = 

1 






index , > 

:.y . 

u,v,m,lx,ly, 

lu, 

lv, lm, tf , gl , g2 . g3 , g4 . o5 . g6 


1, 1. 


« , 

# , i, « , 

* , 

*■ *, i. i. *. i, i. 

JL 

2, # , 

# , 

* , 

#. «, «. 

* . 

*. «. *. *. «. #. *, *. #. 

* 

3, #, 


♦ . 


#, 

*. *. *. «. «, *. #, #. «. 

* 

4 , « , 

* , 

« , 

#, #, #, 

♦ , 


* 


MPMM3D 


hx 

* i 





hy 

* i 





hz 

* 1 





ex 

* i 





*y 

* i 





NODE 

= 3 





SEL 

= 3 





INDEX, X, Y, 

Z,U, V, 

fc'.M 

. LX, LY, L2, 

.LU, 

I. 

i i 

1 *, 

#, 

* 1 •. 

#, 

2, 

* . « , 

*, *, 

• 

« * . * 

* , 

3, 

• • 

V • , 

• . 

« «. t 

1 . 

4, 

f, *, 

« « 

« , 

*, * , (, 

• 


LV.LW, LK.TF, Cl 

G2 

G3 

G4 

G5.G6 

G7 

GF 

G? 

*, *. «. 4. 1, 

1 

1, 

«, 

», • 

• , 

• , 

1 

• «, * * 4, 

* 

«, 

1 , 

*. 1. 

« 

« 

$ . 

«. *. #, t, « 

« , 

• 

#, 

• . 1 

1 , 

« , 

$ , 

*, *, *, «, 1, 

♦ , 

*, 

• 

i, «, 

« , 

• r 

« 


1 1 
* * 
• • 

• • 


1 


♦ 
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MPM2D Iteration 
Info to Screen 


MPM3D Iteration 
Info to Screen 


The important number items are “a,” “ex," and “ey* with “x," “y," 
"m," “tf * “g4,” “g5,” and “g6” on the first line only All other 
numbers are read by the program but not used. The “x” and “y” 
coordinates are used only to compute the true anomaly angle that the 
burn begins at The only mass value remembered is the initial 
mass value. The mass costate is used to scale the constant 
Lagrange multipliers “g4,” “g5," and “g6” in a manner consistent 
with patching the burns together; otherwise, it is not used. 

Listed below is sample screen output from “MPMM2D.F” 


Cur Norm 

It* 

B#st Norm 

Utl « 

Short Tim* 

Bn* 

Bit Writ El 

El* 

0 45051E-00 

1 

0. 45051E-00 

1 

0 . 11269E+G1 

4 

0 3G952E*C0 

1 8 

0 45C51E*00 

45 

C .45051E*00 

15 

0 . 112 69E* 01 

4 

0 . 30952E-0C 

16 

C . 63552E-02 

90 

0 . 63 550E-02 

89 

0 10551E*C1 

4 

0 . 284 ?1E- 02 

1 4 

0 465^SE-02 

135 

0.46575E-C2 

105 

0 . 10606E* 01 

4 

C.24346E-C2 

14 


d • Function EVals * 172 


Total Bum Tima * €.514026424486 
Final Mass * 4.066434637841 
Shortast Bum length * 1 . 1288B3678329 
Shortest Bum is «4 


The first block of text is the iteration table. The column “Cur. 

Norm shows the current 2-norm of the constraint errors in the 
absolute sense. The iteration, or number of times called, at which 
this value was computed is listed in column “It#.“ The lowest norm 
of constraint errors yet computed, next to the iteration number it was 
computed at is given under the “Best Norm (at) #” column. The 
length of the shortest burn at the current iteration is under “Short 
Time“ and the burn with this length is indicated under the “Bn#“ 
column. Finally, the largest absolute value of a constraint 
component for the best norm is listed under “Bst Wrst El.“ with 
“El#“ listing which constraint component this is. 

The iteration table from MPM3D is slightly different. It has the 
following header: 

I ' A7) * TOST c. EL. EL* BST WRST EL EL* 

where “WRST C. EL.” indicates the worst element of the current 
iteration constraint error vector. 

For MPMM2D and MPMM3D, below the iteration table is the 
number of function calls required to reach an error level indicated 
by the tolerance. After this, some statistics of the solution are given. 
The ^Total Bum Time“ is the total amount of time the motor is on. 
The “Final Mass” is the mass of the spacecraft at the end of the 
transfer. The “Shortest Burn Length” is length in time of the 
quickest burn. Finally, the burn number for this quickest, or 
shortest, burn is listed. 
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Data File Output 


IV*2. The Structure of the 
MPMM2D (MPMM3D) Code 


Page 9 


The subroutine MPM2D (MPM3D), if desired, creates an output file 
that gives the status of iterations. The file is named 
‘MPM2D.ISTAT ("MPM3D.ISTAT”). This file is useful for 
computer systems that operate under a queuing system because such 
a system often does not show output to the screen until after execution 
is completed However, such queuing systems usually allow files 
that are created and closed to appear in the users directory. 
Therefore, during execution under a queuing system, the user may 
list the contents of “MPM2D.ISTAT” (“M PM 3 D . I STAT” ) and see' 
current iteration information. The content of “MPM2D.ISTAT” 

( MPM3D.ISTAT”) is three lines long: the first two lines are the 
table headings from the iteration table, the third line is the current 
entry in the iteration table. 

Both the main routine MPMM2D (MPMM3D) and MPM2D 
(MPM3D) contribute to a file named “MPM2D. REPORT" 

(“MPM3D. REPORT”). The first lines in this file gives feedback 
from MPMM2D (MPMM3D) while reading “MPM2D. GUESS” 
(“MPM3D. GUESS”) so that any errors in that file may be easily 
identified. 

The first eleven lines give the header parameters. At the beginning 
of each line, the text from “MPM2D. GUESS” (“MPM3D. GUESS") is 

given, then the number read from that line, and finally, in 

parentheses, the name of the variable which MPMM2D (MPMM3D) 
has assigned this number to. This same pattern is continued as 

MPMM2D (MPMM3D) reads the orbital elements of the transfer 
orbits. 

The twelfth line and lines below are printed as each line of the input 
are read. Following this is a listing of the values of each variable 
used by MPM2D (MPM3D) for the first iteration; then a listing of the 
constraint values when given these variables. 

Next is the iteration table as printed to the screen. Following this a 
total number of calls to MPM2D (MPM3D). Then a listing of 
variables and constraint evaluations for the solution. Finally, at 
the bottom of the file is the solution summary statistics just as 
printed to the screen. 

The other file created by MPMM2D (MPMM3D) is “MPM2D.SOL” 
(“MPM3D.SOL”), the solution file. This file contains the solution to 
the orbit transfer problem in the PAT2D (PAT3D) format. 

The structure of the MPMM2D (MPMM3D) program is generalized 
m the follownng diagram, not intended as a formal flow chart 
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MPMM2D/MPMM3D Diiy nua 



Calls the “dtidimensional nonlinear equation 

r“MPMinrnipef? NF ’ Wlth the guess from “MPM2D GUESS” 

( MPMSD.GUESS ) The solver calls MPM2D (MPM3D) iteratively 

to solve the problem and to numerically compute partial derivatives 
Th, s recurrent use of MPM2D (MPM3D) is illustrated in the 
diagram by a loop with an arrow on it, connecting the two blocks 

MPM2D (MPM3D) evaluates the MPM conditions given the 
variables. For each bum in the orbit transfer problem variables 
* 1 B n URN. This subroutine integrates each burn ic by 
"ST* evaluates boundary conditions for that burn by 

LSOr)F BCC (BC( ; ) ' , [ he derivatives for integration, required for ' 

S’r SUPP ied by FBURN - FBURN is cal,ed repeatedly by 
LSODE during solution of each bum’s initial value problem 
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V. The Patched Method in Two Dimensions (PAT2D) 

The subroutine FUNC is a realization of the Patched Method in two 
dimensions. The file "PAT2D.f' contains an implementation of 
FUNC with the conjugate gradient method. The conjugate gradient 
algorithm was taken from "Numerical Recipes" and is only 
slightly modified from what is presented there. 

PAT2D requires two input files for execution These files specify 
iteration parameters ( "PATCH2D.TOLS" ) and the initial solution 
guess ("PATCH2D. GUESS"). The "PATCH2D GUESS" file must 
be in the PAT2D format (see Appendix A). The format for 
"P ATCH2D.TOLS" is much simpler and demonstrated in the 
example below: 


V.l l Tsing PAT2D to 
Compute Sub- 

Optimal/Extrema] Solutions 


FTOL 

= 

1 . 00000000 0000000 000 00E- 08 

LTOL 

c 

1.000000000000 00 00000 0E- 07 

GTOL 

e 

1 . 00000000 OOOOOOOOOOOOE- 03 

TOL2 

s 

1 . 00 0 00 0 00 OOOOOOOOOOOOE- 05 

ITMX 

= 

200 

MFUN 

= 

200 

Min; 

S 

1000 

ITME 

= 

15 


The FORMAT edit descriptors for the first four lines, containing 
REAL values, are (1X,A6,D27.20) and likewise for the last four 
lines, containing INTEGER values, (1X,A6, 16). The value for 
FTOL specifies the function value stopping criterion, when the 
change in total burn time after a line search is less than FTOL the 
iteration stops. The value for LTOL is the line search tolerance 
GTOL specifies how small the 2-norm of the gradient should be fore 
stopping. TOL2 is the tolerance for DCNLP one-bum solutions 
ITMX is the maximum number of allowed conjugate gradient 
iterations. MFUN limits function calls and MITN limits the 
overall iteration count for DCNLP. ITNB limits the number of 
multiple-shooting iterations performed by BOUNDSCO 

V.2. How PAT2D Works The diagram below shows the general structure of the code in the file 

“PAT2D.f.” 


PAT2D Diagram 
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The Bubroutine FUNC is the heart of PAT2D. This is the function 
that, given the choice of intermediate orbital elements, calculates 
the total bum time for the transfer. FRPRMN is the conjugate 
gradient routine, from “Numerical Recipes,” that iteratively calls 
FUNC and DFUNC (gradient routine) to find the optimal choice of 
intermediate transfer orbits. 


PAT2D has a two-loop structure; there is an inner loop 
(FUNC/ONEBRN) and an outer loop (FRPRMN). The outer loop 
successively changes the transfer orbits until a minimum is found 
in the total burn time (maximum of final mass). The inner loop 
solves the one bum trajectories between each transfer orbit. Solving 
this trajectories yields the bum time s for each intermediate 
transfer. These bum times are summed, giving the output of 

T T XT P * 


Note that each successful outer loop iteration produces a suboptimal 
transfer. This transfer satisfies all the conditions on the state but is 
not an extremal transfer. 


The main routine loads the solution guess and calls FUNC once, 
before FRPRMN does. This is done because there is no assurance 
that the trajectory guesses in the PATCH2D. GUESS file will 
successfully produce a suboptimal solution. The output from this 

“ PATCH 2D INITIAL” and is often a good guess 
for MPMM2D. However, if this is a poor guess, then a good strategy 
is to allow PAT2D several iterations to produce a transfer closer to 
the solution. 


The inner loop iterations are a little complicated. This is the result 
o an attempt to make them robust. It is also designed so that each 
successful inner loop iteration produces a solution to the Two Point 
Boundary Value Problem (TPBVP) with BOUNDSCO, a multiple- 
point shooting algorithm (MS). However, it is widely known that 
direct methods often have a large region of convergence than 
mdirect methods. Therefore, Direct Collocation with Nonlinear 
Programming (DCNLP) has also been implemented. 


The following diagram shows how the ONEBRN subroutine 
interprets the user’s selection as to what is the appropriate first 
action, use MS or DCNLP first? 
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ONEBRN Flow Chart part 1 
(abridged) 


ONEBRN Flo v Chart part 2 
(abridged ) 



Note that a MS guess can be given for DCNLP in this structure A 
DCNLP guess cannot be given for MS because a DCNLP solution is 
required in the conversion process from DCNLP information to MS 
information. 


The next diagrams shows how MS (BNDSCO) and DCNLP (IMSL's 
DNOONF) are incorporated: 



Attempts with either method have a similar structure. If a failure in 
iterations occurs, the guess is perturbed and the method attempted 
again. After each failure, the perturbation size is increased. If MS 
fails too many times, control is handed over to DCNLP. However, if 
DCNLP fails too many times there is no backup and an error exit 
occurs. 


After ONEBRN succeeds in computing a MS solution, the SEL 
parameter is set to 2 for that burn. 
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Output files: 








"PATCH2D HIST" (iteration data) 

PATCH2D. INITIAL (first suboptimal sol) first optimal 
solution obtained, in patch2d format 

"PATCH2D.SOL" contains the extremal solution obtained to 
tolerance 

PATCH2D.BURN (iteration status) prints iteration status, 
file is useful when program is being run under a queuing 
system and screen output is withheld. Printed after a burn is 
solved. 


"PATCH2D.COST" (iteration status); file is useful when 
program is being run under a queuing sys and screen output is 
withheld. Printed after a complete transfer is solved. 
PATCH2D. CURRENT" contains current suboptimal 
trajectory, unless it is the best. 


rAiL,n^u.tstbr 


contains Dest suboptimal trajectory to date 
PATCH2D.PERT" gives information as to the progress of 
solving the current burn. 

FRPRMN. OUT" output from conjugate gradient routine, 

F R P R M N 


"FRPRMN. ITERATES" current output from FRPRMN, for info 
when using a queuing system 
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VL The Multiple Shooting Approach (BND3D) 

The BND3D program implements the modified multiple-point (MS) 
algorithm of BOUNDSCO (Boundary value problem solver with 
Switching Conditions). BOUNDSCO makes use of Newton’s 
method, a Broyden update, and Deuflhard’s relaxation strategy. 

One should refer to the BOUNDSCO manual^ for detailed 
information on BOUNDSCO. Note that BOUNDSCO does not make 
use of an analytical gradient. 

BND3D also has a homotopy loop around BNDSCO. A homotopy 
variable U is defined such that, as the loop repeats, U will change 
from 1 to UMIN (The choice of UMIN is set by the user, but usually is 
chosen as 0). Certain parameters for the orbit transfer problem 
definition are included in the homotopy loop and vary as the value of 
U changes. A tutorial using homotopy is included in the Tutorials 
section. 

The code MP2BND will convert MPMM3D input files into BND3D 
input files. 


VI. 1. Using BND3D to 
Compute Solutions 

BND3D requires two input files: “BND3D. SCRIPT” which contains 
instructions and parameters, and another file (named by user ) 
which contains the solution guess. 

The format of the file “BND3D.SCRIPT” depends on how BND3D is 
to be used. This format is best described line-by-line. The 
character in the first column of each line is ignored. 

The four different layouts of the “BND3D.SCRIPT” file are 
described below’; 


Normal Execution: 
Free Final Time, No 
Homotopy 


• Line 7. (1X,A28) On this line, the name of the file containing the 

solution guess is specified. No more than 28 characters 
are allowed. 

• Line 2. (IX, 16) Here, a “1" indicates that boundary condition 

errors should be displayed to the screen, in addition to the 
normal BNDSCO iteration output; a "0" indicates 
otherwise. Usually, one would place a “0” here; this 
output is usually only useful in finding errors in the input 
file. 

• Line 3: (IX, 16) A “1” on this lines chooses the free final time 

option. 

• Line 4. (IX, 16) A “0” deselects the homotopy option. 

• Line 5. (IX, 16) A 1 on this line tells BNDSCO to insert nodes 

for the switching times in the output; a “0" says not to 


lOberle, H.J. Grimm, W„ “BNDSCO: A Program for the Numerical Solution of Optimal Control 
Problems,” English Translation of DFVLR-Mitt. 85-05. 
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Line 6: (A.D12.5) The value on this line sets the BNDSCO 
parameter FCMIN. FCMIN is the lower limit of the 
relaxation factor. 

Line 7: (A,D12.5) The value on this line sets the BNDSCO 
iteration tolerance. 

Line 8, (IX, 14) The maximum number of iterations. 

Line 9: (1X.A28) The name for the file containing the solution 

Line 10: (IX, 16) A “1” on this line requests detailed solution 

information (“BND3D. EXTRA” and the file named on 
the next line). A “0” indicates otherwise. 

Line 11 ^ (1X.A28) The file name for additional information (if a 
“1” on the previous line). 


Fixed Final Time; 
No Homotopy 


• Line 1: (1X.A28) On this line, the name of the file containing the 
solution guess is specified. No more than 28 characters 
are allowed. 

Line 2. (IX, 16) Here, a 1 indicates that boundary condition 

errors should be displayed to the screen, in addition to the 
normal BNDSCO iteration output; a "0” indicates 
otherwise. Usually, one would place a “0” here; this 

output is usually only useful in finding errors in the input 
file. 
















Line 3: (IX, 16) A “0” on this lines chooses the fixed final time 
option. 

Line 4: (A.D12.5) The value for the final time. 

Line 5: (IX, 16) A “0” deselects the homotopy option. 

Line 6: (IX, 16) A “1” on this line tells BNDSCO to insert nodes 
for the switching times in the output; a “0” says not to. 

Line 7: (A.D12.5) The value on this line sets the BNDSCO 
parameter FCMIN. FCMIN is the lower limit of the 
relaxation factor. 

Line 8: (A.D12 5) The value on this line sets the BNDSCO 
iteration tolerance. 

Line 9: (IX, 14) The maximum number of iterations. 

Line 10: (1X,A28) The name for the file containing the solution 

Line 11: (IX, 16) A “1” on this line requests detailed solution 

information (“BND3D. EXTRA” and the file named on 
the next line). A M 0 W indicates otherwise. 

Line 12 . (1X.A28) The file name for additional information (if a 
1 on the previous line). 


Free Final Time, 
Homotopy Activated 


• Line 1: (1X.A28) On this line, the name of the file containing the 

solution guess is specified. No more than 28 characters 
are allowed. 

• Line 2: (1X.I6) Here, a “1” indicates that boundary condition 

errors should be displayed to the screen, in addition to the 
normal BNDSCO iteration output; a “0” indicates 
otherwise. Usually, one would place a “0” here' this 

output is usually only useful in finding errors in the input 
file. r 

Line 3: (IX, 16) A “1” on this lines chooses the free final time 
option. 
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• Line 4: ( IX, 16) A “1" selects the homotopy option 

• Line 5: (IX, 16) the suggested number of homotopy loops to 

perform 

• Line 6 (*) Enter UMIN, the value of the homotopy variable to stop 

at. The homotopy variable, U, starts at 1 and ends at 
UMIN. Enter “0.0” here to attempt to achieve the values 
below. 

• Line 7: (*) Enter the desired maximum thrust level 

• Line 8. (*) Enter the desired specific impulse 

• Line 9 (*) Enter the desired final orbit semimajor axis 

• Line 10: (*) Enter the desired final orbit eccentricity 

• Line 11: (*) Enter the desired final orbit argument of perigee 

• Line 12: (*) Enter the desired initial orbit semimajor axis 

• Line 13: (*) Enter the desired initial orbit eccentricity 

• Line 14: (*) Enter the desired initial orbit argument of perigee 

• Line 15 <*) Enter the desired initial orbit argument inclination 

• Line 16: ( IX, 16) A “1” on this line tells BNDSCO to insert nodes 

for the switching times in the output; a “0” says not to 

• Line 17: (A,D12.5) The value on this line sets the BNDSCO 

parameter FCMIN. FCMIN is the lower limit of the 
relaxation factor. 

• Line 18 : (A,D12.5) The value on this line sets the BNDSCO 

iteration tolerance. 

• Line 19: (1X,I4) The maximum number of iterations. 

• Line 20: (lX r A28) The name for the file containing the solution 

• Line 21: (IX, 16) A “1” on this line requests detailed solution 

information (“BND3D. EXTRA” and the file named on 
the next line). A “0” indicates otherwise. 

• Line 22: (1X,A28) The file name for additional information (if a 

“1” on the previous line). 


• Line 1: (1X,A28) On this line, the name of the file containing the 

solution guess is specified. No more than 28 characters 
are allowed. 

• Line 2: (IX, 16) Here, a “1” indicates that boundary condition 

errors should be displayed to the screen, in addition to the 
normal BNDSCO iteration output; a “0” indicates 
otherwise. Usually, one would place a “0” here; this 
output is usually only useful in finding errors in the input 
file. 

• Line 3: (IX, 16) A “0” on this lines chooses the fixed final time 

option. 

• Line 4: (A,D12.5) The value for the final time. 

• Line 5: (IX, 16) A “1” selects the homotopy option. 

• Line 6: (IX, 16) the suggested number of homotopy loops to 

perform 

• Line 7: (*) Enter UMIN, the value of the homotopy variable to stop 

at. The homotopy variable, U, starts at 1 and ends at 
UMIN. Enter “0.0” here to attempt to achieve the values 
below. 

• Line 8: (*) Enter the desired maximum thrust level 

• Line 9: (*) Enter the desired specific impulse 

• Line 10: (*) Enter the desired final orbit semimajor axis 


Fixed Final Time, 
Homotopy Activated 

(in this case, the 
fixed final time is 
also achieved 
through the homotopy 
loop) 
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Line 11 : 
Line 12. 
Line 13: 
Line 14: 
Line 15: 
Line 16: 
Line 17: 


(*) 

(*) 

(*) 

(*) 

(*) 

(*) 


Enter the desired final orbit eccentricity 
Enter the desired final orbit argument of perigee 
Enter the desired initial orbit semimajor axis 
Enter the desired initial orbit eccentricity 
Enter the desired initial orbit argument of perigee 

(IX d !t irC ?- initial ° rbit ar ^ ment inclination 

<1X,I6) A 1 on this line tells BNDSCO to insert nodes 

Line 1R anm ^ in the ° utput; a “°" sa ys not to 
8. (A.D12.5) The value on this line 6ets the BNDSCO 

parameter FCMIN. FCMIN is the lower limit of the 
relaxation factor. 

Line 19: (A.D12.5) The value on this line sets the BNDSCO 
iteration tolerance. 

I?; ma ” mom ”“">>>« of iterations. 

lZ 22 "T r r ,he containing the solut.on 

. (1X46) A 1 on this line requests detailed solution 
information (“BNDaD.EXTRA" and the file named on 
r „ he ™ Xthne) ^ 0 indicates otherwise. 
me 23: (DU28) The file name for additional information (if a 
l on the previous line). 


format The fin t T” ,n “BND3D.SCRIPT") has a specific 

tranSfCr P aramet ers of 

(S!a9 F30 if) 5^ 8 h8Ve F0RMAT ^it descriptors 

(1X.A9.F30.15). These parameters are as follows and in this order 


MU 

REQ 

J2 


GO 

BETA 

RO 

ROU 

S 

CD 

ISP 

THRUST 

AI 

El 

OMEGAI 

RAI 

M 

AF 

EF 

OMEGAF 

RAF 

I-F 


gravitational constant of the central body (1 0 for no 
dimensions) 


equatorial radius of the central body 
constant describing the mass distribution of the 
central body; for Earth t/ 2 =1082.61xl0*® 
acceleration at sea-level 


constant from the atmosphere model describing air 
density variation in the prescribed altitude region 

f q, ^ 


atmosphere density at the altitude 
cross-sectional area of the craft 
drag coefficient 
specific impulse 
maximum thrust 
initial semimajor axis 
initial eccentricity 
initial argument of perigee (degrees) 
initial right ascension (degrees) 
initial inclination (degrees) 
final semimajor axis 
final eccentricity 

final argument of perigee (degrees) 
final right ascension (degrees) 
final inclination (degrees) 
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The 21st line (IX, 15) gives the number of intervals (# nodes • 1). 

The next line is a dummy string line (1X,A) that, on output, is used 
to provide a header for the data in the following lines (useful in 
plotting results). 

The next (# nodes) lines gives the BND3D state at each node with 
edit descriptors (1X,F30.15,25(A2,F30.15)). The BND3D state is as 
follows 

0 1 2 3 < 5 6 7 e 9 10 11 12 13 14 15 

(T. X, Y. Z, U. V. W, M , L-X, L-Y, L-Z. L-U. L-V. L-K, L-K, TF 

( FINAL ORBIT ) { INITIAL ORBIT } 

16 17 18 19 20 21 22 23 24 25 
Gl, G2 , G3 , G4, G5, G6 , G7 . G8. G9 , G10) 

<X,Y.Z> IS POSITION <L-X.L-y.L-Z> IS LAKBDA-P 

<U , V , W> IS VELOCITY <L-U, L-V. L-W> IS LAMBDA -V 

M IS MASS 

L-M IS LAMEDA-K 

T IS THE NORMALIZED TIKE [0.1] 

Where TF is the final time and G# are components of the constant 
Lagrange multipliers (v); G1-G5 being v for the final boundary 
conditions and G6-G10 being v for the initial boundary’ conditions 

The nodes are entered in the reverse order, starting with the final 
node and ending with the initial node. 

Following the node information is a line (IX, 15) for the number of 
switching points. It is suggested to use an even number of switching 
points - this indicates to BNDSCO that the first and last intervals are 
burn arcs 

The next lines (1X,F30.15), one for each switching point, give the 
switching times in normalized time [0,1]. No lines after these are 
read. 


VL3. How BND3D Works BND3D supplies the necessary routines (F and CON) to BNDSCO 

“F” supplies the derivatives of the state and “CON” evaluates the 
boundary conditions. The routine “BCC" computes repeated 
formulas, “LSG” loads the solution guess, “SAVSOL” saves solution 
data in the same format as the guess data. The routine “DIFSYB" 
performs numerical integration. 
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d i agram be,ow indicates the interdependence 
BND3D subroutines. 


of the 
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VII. The Minimizing Boundary Condition Method (MBCM3D) 

The Minimizing Boundary Condition Method (MBCM) is a relaxed 
simple shooting algorithm. Instead of using a multidimensional 
nonlinear equation solver for the two point boundary value problem 
(TPBVP), it transforms the TPBVP into a nonlinear programming 
(NLP) problem. 

As included in ORBPACK, MBCM3D uses the square of the 
Hamiltonian as the NLP cost function. All other boundary 
conditions are taken as NLP constraints. 


VII.l. Using MBCM3D to 
Compute Solutions 


MBCM3D requires one input file, MBCM3D. GUESS This file has a 
very specific format. The first 47 lines of this file have the 
FORMAT edit descriptors (1X,A9,E30.15). They describe, in the 
following order: 

MU gravitational constant of the centra] body (1.0 for no 

dimensions) 

REQ equatorial radius of the central body 

«J2 constant describing the mass distribution of the 

central body; for Earth «7 2 =1082.61xl0' 6 
GO acceleration at 6ea-level 

BETA constant from the atmosphere model describing air 

density variation in the prescribed altitude region 
RO Tfy +REQ 

ROU atmosphere density at the altitude r & 

S cross-sectional area of the craft 

CD drag coefficient 

ISP specific impulse 

THRUST maximum thrust 

AI initial semimajor axis 

El initial eccentricity 

OMEGAI initial argument of perigee (degrees) 

RA1 initial right ascension (degrees) 

I-I initial inclination (degrees) 

AF final semimajor axis 

EF final eccentricity 

OMEGAF final argument of perigee (degrees) 

RAF final right ascension (degrees) 

I-F final inclination (degrees) 

[the next 14 lines give the initial state] 

TF transfer time 

[the next 10 lines give G1-G10] 

ACC solution tolerance 


Where G# are components of the constant Lagrange multipliers 
(v); G1-G5 being v for the final boundary conditions and G6-G10 
being v for the initial boundary conditions. 

The last line of “MBCM3D. GUESS” (1X,A9,I10) gives the maximum 
number of iterations. 
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VH. 2 . How MBCM3D Works 


MB CM3 D Flow Diag ram 


SSnS-'®'::! * BND3D file named 

BND3D. GUESS into a MBCM3D guess file ( tt MBCM3D GUESS"). 

MBCM3D uses VF02AD to solve the NLP problem. VF02AD uses 
reverse communication: the main routine calls OF to compute NLP 
cost and c °™traints given input; then GRD to compute gradients 
then calls VF02AD to compute the new iterates. The main routine 

VF02A S n heSe , neW lt€ratCS aS input for 0F and «Peats the loop until 
VF02AD signals convergence. 

OF evaluates the TPBVP as a NLP. The shooting problem is 

TZZtr th R w a Runge * Kutta integration routine. Integration 
of the shooting problem is interrupted often to check the sign ofthe 

„i™? 8 TTa ’ f a ,' ig " Ch “" “ «•» integration 

interval is adjusted until the exact switching point is located 

Dunng this process, OF keeps track ofthe sign ofthe switching 
T “ d a PP r °Pnately adheres to the optimal switching law. 
This should ensure that the switching law is followed, however it is 

always prudent to check the switching law after a solution is 
claimed. 

MBCM»X a ontto« W ‘ ndiCa " S ‘ he * nter< i e P en( i ence of the 


MAIN 




o 

^3 


V 




Dj 


n VKJ2AD ] 


IBCC 


INTEGRATION 

LOOP 
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VJLil. Tutorials 


Vn.l. Planar Five Burn 
Transfer 


Use GSHOOT to 
Construct a Guess 


The following tutorials demonstrate some aspects of using 
ORBPACK that the user may commonly encounter. 

This tutorial demonstrates the use of the supplied code in solving a 
planar transfer from a circular LEO to circular GEO. The initial 
radius is 6600 km, the final radius is 42241 km. The initial rocket 
motor thrust is 9.918 kN; its I g p is 450 seconds. The initial mass is 
20980 kg. A five burn solution is desired. 

After nondimensionalization, these parameters are: initial 
mass=10, thrust=0.5166, go=l, Isp=0.5673, initial radius=l, final 
radius=6.4 

Based on the characteristics of these types of transfers, the following 
guess for the transfer orbits may 6eem reasonable 


a 

e 

1.285 

0.2189 

1.570 

0.3584 

i 1.856 

0.4550 

3.707 

0.7262 


All their apses are aligned and the final transfer orbit is similar to 
the Hohmann transfer orbit. 

The trajectory for each bum will now be guessed using GSHOOT. 
The "INDIRECT.DAT” files produced by GSHOOT will then be 
concatenated together to form an “MPM2D. guess” file The first 
burn input file for GSHOOT (“GINPUT”) is supplied as 
“Tutorials/2D 5burn/GSHOOT/burn 1/GINPUT” and listed below: 


Mu 

» 

1 .00 

Go 

£ 

1 .00 

Isp 

£ 

C . 567J 

Thrust 

£ 

0.5166 

Me 

X 

10.0000 

ec 

£ 

I .00000 

e: 

X 

c.oc: 

wc 


o.ooc 

ad 


1.265 

ed 

£ 

0.219 

wd 

£ 

0.000 

TNAX 

* 

O.OOC 

NGS 

£ 

10C 

NIX 

£ 

3 


GSHOOT reports: 


constant Lagrange multipliers (initial l 
C 0 15245E* CO 0.96620E-0C 0.14561E-C1 

Best initial true anomaly 
vc* 0.53047E-C: 

Best transfer time 
tf* 0 . 19312E+01 

Best relative errors (h,ex,ey,Hs) 

! G .. 0.1B817E-0P - 0 . 49820F-C1 0.15555E-02 0.275P9E-02 

The resulting file has been supplied as “Tutorial s/2 D 
5burn/GSHOOT/burn 1/INDIRECT. DAT" The second burn 
“GINPUT” is rTutorials/2D 5burn/GSHOOT/burn 2/"]: 
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Mu 

- 1.00 

Go 

■ 1.00 

lsp 

- 0.5673 

Thrust 

« 0.5166 

Mo 

* 10.0000 

ao 

* 12 05 

ec 

• 0.219 

wo 

■ 0.000 

ad 

« 1.570 

ed 

- 0.3584 

wd 

■ 0.000 

TMAX 

- 0 COO 

NGS 

- 10 C 

NIX 

■ 3 


GSHOOT reports: 


Best constant Lagrange multipliers (initial) 

C. . 0 7C359E*00 C.17901E-00 -0.244C2E-14 

Best initial true anonaly 
vo- 0 . 56451E*01 
Best transfer time 
tf- 0 . 1145BE-01 

Best relative errors (h,ex,ey,Hs) 

C.. 0 . 10846E-07 -0.82645E-02 -0.16307E-02 0.32135E-03 

The resulting file has been supplied as “Tutorials/2D 
5burn/GSHOOT/burn 2/INDIRECT.DAT“ The third burn •" is 
(“Tutorials/2D 5burn/GSHOOT/burn 3 /“)• 


Mu 

■ 1.0C 

Go 

■1.00 

lip 

■ 0.5673 

Thrust 

* 0.5166 

Me 

* 10.0000 

aO 

* 1.570 

eO 

* 0.3504 

w0 

* 0-0CQ 

ad 

* 1.856 

ed 

■ 0.4550 

wd 

* 0.000 

TMAX 

* 0.000 

NGS 

* 100 

NIX 

* ? 


GSHOOT reports: 


Best constant Lagrange multipliers (initial) 
c . .. 0.54451E*00 0.26192E*C0 -0 . 1033CE-14 

Best initial true anomaly 
vo- 0.6C064E-01 
Best transfer time 
tf« 0 . 79429E+00 

Best relative errors (h,ex,ey,Hs) 

G. . 0 . 925“ 7 4E-08 C.48454E-02 0.13266E-01 - 0 . 35436E- 02 


The resulting file has been supplied as “Tutorial s/2D 
5burn/GSHOOT/burn3/“ The fourth burn (“Tutorials/2D 
5burnyGSHOOT/burn 4/ u j: 


mu 


1.00 

Go 


1.00 

lsp 


0.5673 

Thrust 


0.5166 

Mo 

- 

10.0000 

ac 


1 656 

eo 


0.4550 

wo 

- 

0.000 

ad 


3 707 

•d 


0.7262 

wd 


0 .000 

TMAX 


0.000 

NGS 


100 

NIX 


3 


GSHOOT reports: 
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Best constant Lagrange multipliers (initial) 
c 0 . 44412E* 00 0.3C915E-00 0.35928E-14 

Best initial true anomaly 
vo* 0.53^82E-01 
Best transfer time 
tf* 0 . 1 826 5E*0l 

Best relative errors (h,ex,ey,Hs} 

G 0 56826E- OF -0.39904E-03 0 17988E-02 -0.36833E-C2 

The resulting file has been supplied as w Tutorials/2D 
5burn/GSHOOT/burn AT The fifth burn [“Tutonals/2D 
5burn/GSHOOT/burn 5/"): 


mu * ].0C 

Go *1.00 

2sp * 0.5673 

Thrust » 0.5166 

Mo * 10.CC00 

ao * 3.7C7 

* 0.7262 
wc *0.000 

ad =6 400 

ed & 0.0CC0 

wd * C . 0 0 0 

TKAX * 0.000 

NGS * 100 

MX * 3 


GSHOOT reports: 


Best constant Lagrange multipliers (initial) 

C. - 0 . 28C15£*0C -0.718 C2E*00 -0.€37;5E*0C 

Best initial true anonaly 
vo* 0 . 30096E- 03 
Best transfer time 
tf= 0 . 3221 9E* 01 

Best relative errors (h.«x,ey,Hs) 

G... 0 . 26C77E- ii -0 . 93204E-02 -0.25981E-01 0.53e-DB£-01 

The GSHOOT output has been supplied as "Tutorials/2D 
5burn/GSHOOT/burn 5/“ 


The files easily concatenate. The resulting file has been supplied 
as “Tutorials/2D 5burn/GSHOOT/MPM2D.guess“ 


Attempt Computation 
of Solution with 
MPMM2D 


At this point, we have a solution guess for the entire trajectory’ in the 
PATCH2D format. One option for obtaining the solution is to run 

MPMM2D with this input. However, one may get a iteration history 
like this: 


MPMM2D CWnnt 
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Uf B * ,t N0nr Ut) • S*°rt Tim# Bnf Bsc Wrst El 


0.«S735E*01 i 0 68735E-01 1 

0 6eT35E*0l 45 0.68735E*01 45 

0 6B73SE-01 90 0 6B735E-01 45 

BuC Fossibl* conflict in orbit choice 
A*-2 617712643152 
E»2 .335666952254 
W«2 556150238017 
(LOCATION «l] 

BCC. Possible conflict in orbit choice 
A«-2 617712643152 
E*2 335666952254 

W*2 .556150238017 
{LOCATION fl] 

B'JRN WARNING. BCC CLAIMS AN ERROP 

IN THE INITIAL PCINT CALCULATION 
WI«5 684341 886 06 08E- 14 
W2-1 858576979153 
W3=0 7387094236308 
{LOCATION «1 j 

BCC. Possible conflict in orbit choice 
A*4 117497625609 
E*1 458915419969 
W=-0. 5075614176646 
(LOCATION *1) 

INCONSISTENT 
A* (l#0‘E**2i . LT , 0E0 
STOP (called by BCC > 

CP. 20.155s, Wallclock: 29 935s, 

HWM mer.: 213617, kwm stack. 26610 


0 79429E* 00 
0 79429E* 00 
0 79429E* 00 


3 C.34C45E*:; 
3 0 34C45E-01 
3 0 34C45E-C1 


33 7% of 2 -CPU Machine 
Stack overflows; 0 


Ela 

4} 

43 

43 


Note that the current norm error started at 6.3735; though such a 
large error does not always induce failure of MPMM2D, it may. 


If MPMM 2D Fails, 
Use PATCH2D 


In such a situation, the more robust PATCH2D is useful. Since the 
file format is identical, this is very convenient. PATCH2D does 
require one i additional input file, for its inner loop tolerances. The 
file is called PATCH2D.tols” and for this tutorial, it has been 
supplied as “Tutorials/2D 5burn/PATCH2D/PATCH2D.tols“ and 

listed below-: 


FTOL * 
LTCL * 
GTCL * 
TOL2 * 


1 OOOOCOCOCCCOOCOCOOCOE-C8 
1 . OOOOOCOOCCOCOCOOOOOOE-07 
1 . COOCCOCOCCOOOO:jCOOOE-03 
1 . OOOOOOOOOOOOOOCOOOOOE-05 


We have chosen a rather strict tolerance for “function 
improvement” convergence, a slightly less strict tolerance for “line 
search convergence, a very loose tolerance for “gradient norm" 
convergence, and a rather loose convergence tolerance for DCNLP 
iterations. 


It needs to be said that the drawback to PATCH2D is its speed. For 
this tutorial, PATCH2D was run. After renaming 

MPM2D . GUESS” to “PATCH2D. GUESS” and running PATCH2D, 
we 6ee the following iterations; 
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rune ti or TCL (FT0L> * l.E-F 
Gradient. TOL IGTOL * « 1 E-3 
Line Search TOL ILTOLl « 3 E-7 
Hax * iterates UTMAXi * 2CC 


IT* 

Cost rune Improvement 

Gradient 

Cr ; ter; or. 

Line 

• eva 1 

0 

0 664 55E* 01 O.OOOOGE*OQ 

0 . 4721 IE* C2 




I 

0 664 55E*03 0 OOOOOE-OO 

0 . 4721 IE* C2 

0 66837E-C1 

7 

12 

2 

0 . 66C21E* 03 -0 43438E-01 

0.1443 9E* 02 

0.22 906E- 02 

6 

22 

3 

0 . 6591 IE* 01 -0 54372E-01 

0 . 56984E*C1 

C 52307E-C2 

4 


4 

0 . 65885E*03 -0 56986E-03 

0 . 121 30E*02 

0 1 074 0E- 02 

5 

41 

5 

0 . €58 3 IE* 01 -0.62356E-01 

0 . 622 5CE*Ci 

0 . 55562E-02 

4 

t *. 

6 

0 6S554E*01 -0 90139E-CI 

0 . 1 8560E* 02 

0 13395E-C1 

5 

63 

7 

0.65467E-C1 -0 96636E'C1 

0 . 119S2E* 02 

0 . 12299E-01 

4 


B 

0 . 65425E* 01 -C.1C299E-00 

0 . 39034E* 01 

0. 17B44E-02 

A 

82 

9 

0.6541 6 E* 01 -0.1036 BE- DC 

0 . 75207E* 01 

0 8P954E- C2 

5 

52 

10 

0.65372E*01 -G.10B33E*00 

0 . 71069E* 01 

0 852 66E- 02 

3 

ic: 

11 

0 . 65329E*01 -0.112S9E*00 

0 . 8O960E* 01 

C . 2134 IE- 02 



12 

0 6531BE*01 “0 . 1 13 66 E* 00 

0 . 32664E* 01 

0 . 13259E-C2 

6 


13 

0 . 65312E* 01 -0 . 1 1432E*00 

0 . 13498E* 02 

C. 5005GE-04 

€ 

131 

2 4 

0 . 6531 IE- 01 “C 11435E*00 

0 . 55B07E* 00 

C 31056E-G4 

5 

14C 

15 

0 65311E*C1 “0 . 1 1436E* 00 

0 . 33032E* 00 

0. 1 64 64E- 04 

3 

145 

If 

C . 653 1 IE* 01 -C . 11437E-00 

0 . 76877e*00 

0. 18622E-C3 

3 

1 55 

17 

0 653 1 0E*C1 *0 11 44 6E* 00 

0.11374E*C2 

0 24295E-C3 

3 

16 5 

1 ? 

C . €53C?E*0i *0 . 11456E*0C 

0.19365E*C: 

C . 14 95 3E* 03 

4 


IS 

0.6530eE*Cl -0.11466E-00 

0. B3243E- 00 

G. 1C 03 EE- 03 

4 

18 6 

20 

0 6530BE-C1 -0.11471E*00 

0 . 54440E* 00 

0 . 1 I25BE- 04 

4 

195 

2 : 

0 6 53C6E* 01 -O.i:47lE*O0 

O.43573E*0C 

0.76724E-C4 

t 

206 

22 

0 6 53 C7E-C1 - C . 1 1 4 7SE* 00 

0.80416E*00 

0 1O644E-03 

■j 

215 

23 

0.653 07E* 01 -0 . 114 £1E* DC 

0 . 25569E*02 

0 . 13474E-C3 

4 

2 2 4 

24 

C 6 53 06E* 01 “0. 11467E-00 

0 . 43370E* 0C 

C. 33C22E-04 

3 

233 

25 

C 653 C6E* Cl -0 . 11489E*00 

0 . 70424E* 00 

0 . 24 93 0E- 04 

4 

^ ^ ^ 

2 f 

0.65306E-01 -0 . 214 90E*00 

0.56308E*DC 

0.209S7E-C3 

4 

2 53 

27 

0 653C5E* Cl -C.11501E*00 

0 . 14167E-C1 

C . 272C5E- 03 

4 

262 

2 £ 

0.653 03E* 02 -0.115141*30 

0 . 30533E* 01 

C 10666E-C1 

6 


25 

0.652501*01 *0 . 12046E*00 

0 . 759“1E*C1 

0. B3111E-G2 

7 

285 

30 

0 . 65209E*01 -0 . 12463E* 00 

0 . 21314E*01 

0.34787E-G3 

5 

2oq 

31 

0.652 C7E* 01 *0 . 12481E* 00 

0. 18041E*01 

C . 23 51 5E- 03 

6 

3C9 


The PATCH2D code had been left to run overnight, about 12 hrs. It 
did not satisfy any convergence criterion by the 31st iteration, 
execution was terminated. The output file “PATCH2D.BEST" has 
been put into in the “Tutorial” folder as “Tutorials/2D 
5burn/PATCH2D/PATCH2D.BEST“ 


Use PATCH2D 
Output for MPMM2D 


Now, this file was renamed to “MPM2D. GUESS” and used for input 
to “MPMM2D.” The iterations are listed below: 


CUE. . NORM 

IT* 

BEST NORM 

t AT) # 

SHORT TIME 

BN* 

BS7 WF.S7 EL 

EL* 

0. <C24CE*00 

1 

0 .4C240E*CC 

1 

0.67665E*:: 


C v ■' 


C . 4 024CE* 00 

45 

0 .40240E-00 

41 

0 . 67665E* CO 

3 

C 31525F* CC 

34 

■5 (■ 

0 35362E-C2 

90 

0 . 35361E-C2 

72 

0 . 10959E* Cl 

3 

C 19245E-02 

0 . 24412E-06 

135 

0.3588CE-10 

123 

o.ii2E8E*c: 

4 

C ■sf tc pr _ - - 


0 . 654 24 E- 07 

1 60 

0 .3G394E-10 

174 

0.1126EE*:: 

4 

C 2 4 ~ ” * r - ‘ ' 

I ‘ 

C .29068E-10 

225 

0 . 29C6SE- 1 0 

225 

o . H26EE* ; : 

4 


^ ' 

C.72452E-07 

270 

C .29068E-10 

225 

c . H28BE*:: 

4 

f C*2 0€ r --' 


0 . 11927E-06 

315 

0 . 26966E* 10 

276 

0 . 212E6E-0I 

4 

C 2* * C 


0 .47056E-07 

360 

0 . 22244E-10 

331 

0 . 112 SEE* Cl 

4 

C. 16456E-1C 

26 

0 . 28231E-08 

405 

0 .2D320E-10 

360 

0.U266E-C1 

4 

0.1522 0E- 1 C 

0.23782E-06 

450 

0 . 15318E-10 

441 

0 . 11266E* 01 

4 

C . 1 06 5CE- 1 C 

2 * 

0 . 1 4141E-1 0 

495 

0 . 13615E- 10 

493 

0.U28BE*C1 

4 

0 . 597 6CE- 1 3 
C . 99" 6CE- 11 

2 f 

0 . 11765E-06 

540 

0 . 13615E- 10 

493 

0 . 11266E*01 

4 

2 ( 

0.40513E-06 

565 

0 . 13339E-10 

544 

0 . 11286E*C1 

4 

C 92B66E-1I 

26 

0 . 47061E-07 

630 

0 . 11 624E-10 

599 

0.1126eE*Cl 

4 

C .76525E-11 

0.78723E-07 

675 

0 . 96522E-11 

654 

0.11266E*0: 

4 

c €C75 :e*;i 

26 

*• FATAL 

ERROR 3 from NECWF 
The user may try a 

The iteration has 
new initial guess 

not made good progress 


Obviously, the solution was found; however, a shortcoming in the 
NEQNF solver did not allow it to claim convergence. This seems to 
be common among nonlinear equation solvers. An easy fix is to 
perturb the guess slightly. In this case, the eccentricity of the first 
transfer orbit was perturbed from 

I «* ■ C.14<237S3690e?2B3e260E->CC 
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to 

| *X * 0 16433 ^5 3 6 90672 83 626 0E*O0 


For this new guess, in the “Tutorial” folder as “Tutorials/2D 
5burn/MPM2D. GUESS," the MPMM2D iterations are: 


CUR NORM 

IT* 

BEST NORM 

(ATI * 

SHORT TIME 

BN* 

BST WRST EL. 

EL* 

0 40418E-00 

1 

0 . 4041BE-00 

1 

0 . 67665E-00 

3 

0 . 31525E* 00 

34 

C . 404 1 8E*Q0 

45 

0.40418E*00 

41 

0 67665E-C0 

3 

0 . 31525E-0C 

34 

0 30687E-C1 

90 

0.30687E-01 

56 

0. 10688E-01 

3 

0 . 14737E-01 

26 

0 4 66 3 0E- 07 

135 

0 . 21092E-10 

122 

0 . 1 1268E* 01 

4 

0 . 13635E-1Q 

2 < 

0 6 09C6E- 07 

1 B 0 

0 . 18042E- 10 

172 

0 . 11288E*0i 

4 

C . 14477E-10 

22 

0 . 302 I 4E-06 

225 

0.17B36E-10 

220 

0 . 11288E*01 

4 

0 . 14065E- 1C 

22 


REQUIRED f FUN CT I ON EVA LS * 268 


TOTAL BURN TIKE * 6 . 51 3750674051 
FINAL HASS * 4 . 06 83 8*7 6C501 5 
SHORTEST BURN LENGTH * 1 . 128B31615888 
SHORTEST BURN IS 44 


SOLUTION SAVET 


The solution file is given in the “Tutorial” folder as “Tutorials/2D 
5burn/MPM2D.SOL”. 


VIIJ2. Convert MPMM3D 
File to BND3D File, Run 
BND3D 


This tutorial demonstrates how to use MP2BND to convert a 
MPMM3D file to a BND3D file. 


The file “Tutorials/MPM to BND3D/MPM3D. GUESS" is a solution 
to an orbit transfer problem, as claimed by MPMM3D. The 
particular problem it solves is not relevant, but it will be clarified 
anyway. The header of this file follows: 

TOL = 0 . 10000000000000000000E-08 
MU = 0 . 10000000000000000000E+01 
T = 0 . 51658300000000068053E+00 
Go * 0 . ICOOOOOOOOOOOOOOOOOOE+Ol 
Isp = 0 . 56730999999999909278 E+ 00 
hxo = 0 . 477I5876030000003993E+00 
hyo = O.OOOOOOOOOOOOOOOOOOOOE+OO 
hzo = 0 . 87881711269999840397E+00 
exo « O.OOOOOOOOOOOOOOOOOOOOE+OO 
eyo = O.OOOOOOOOOOOOOOOOOOOOE+OO 
hxf = O.OOOOOOOOOOOOOOOOOOOOE+OO 
hyf = O.OOOOOOOOOOOOOOOOOOOOE+OO 
hzf = 0.25298517739999937248E+01 
exf « O.OOOOOOOOOOOOOOOOOOOOE+OO 
eyf * O.OOOOOOOOOOOOOOOOOOOOE+OO 
NORB * 5 


Run MP2BND 


The orbit transfer is, therefore, from LEO to GEO and circle to circle 
in 6 burns. Now, suppose we want to further investigate this problem 
with the more general BND3D code, so that oblateness and drag 
effects can be modeled. 

The main task here is to simply run MP2BND. This code will 
create the file “BND3D. GUESS” which hds been supplied as 
“Tutorials/MPM to BND3D/BND3D. GUESS.” 
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Run BND3D to check 


It is prudent at this point to use BND3D to check MPMM3D’s results 
In this tutorial, the following "BND3D. SCRIPT" file was used 

BND3D . GUESS 
0 
1 
0 
0 

Id-4 
Id- 1 0 
100 

BKD3D. SOL 
1 

BND3D . REINT 


which is supplied as “Tutorials/MPM to BND3D/BND3D SCRIPT - 
This says the input file is “BND3D. GUESS,” don’t show B.C. errors 
to the screen, solve with free final time, don’t include switching 

P° intS - nodes ln the o^put, FCMIN=lD-4, TOL-lD-10, use no more 
tnan 100 iterations, save solution as “BND3D.SOL,” provide 

“? f ° and 6ave this info in “BND3D.REINT.” The outDut 
BTsD3D produces to the screen is listed below: 

B.C .S ’ C 

FPXI FINAL TIKE I 
HOKCTOPY . 0 

1 . OOCCOODOOOOOOCGCO 
0 . OOOOOOOOOOOOOOOOOOE' 

0 OOOOOOOOOOOOOCOCOOE 

1 . ooooooooooooooooo 

0 . 0000CC00 00000 DC 0C0E 
O.OOOOOOCOOCOCCOCOOOE 
0 . OCOCOOCOOOCCOCOCOCE 
0 . OOOOOOOOOOOOOCOCOOE' 

0 OOOOCOOOCOOOOOOOOOE' 

0.5673C5999999996962 
C .516563000000001014 

1 .OCC0CCC0010536792 
0 . OOCOCOOOOOOCOOOOOOE-OC 
0. OOOC I OCOOOCOOCCCCGGE*OC 
89 . 9999999997066680 
26 5000C00009010ei9 
6 .40014999641091026 

o . occooc::-:cooccooooe*go 
o . oocooco:ooccoooogoe*gc 
c. co:cccooc:o:occoooe-c: 
o . oocoooooooooocooooe^oo 

NOTE: AN 3LES MUST BE IN DEGREES 
K= <4 
*N= 2 5 


KU- 
REw* 
J2 * 

GO* 
BETA= 
RO= 
t yj- 
s= 

CD* 
ISP* 
THRUST 
AI = 

El* 

OMEGA I 
KAI = 
1-1 = 
AJF= 

EF* 

OKEGAF 

RAF* 

2 -F = 


*00 

*00 

00 

CO 

CO 

00 

oc 


INITIAL DATA 


PRESCRIBES REU.TIVE PRECISION .10D-0» 
MAXIKUK PERMITTED NUMBER OF ITERATIOKSIOC 


AES ERR 


LEVEL 1 


LEVEL2 


LT.TL3 RELAX NEW CON’D (Mi 


N0FJ-: 
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Useful Information 
in BND3D. EXTRA 


0 

. 14D-07 

. 1 ID- 07 


.14D-07 

. 1 ID- 07 

1 

. 14D-07 

.67D-07 


. 14D-07 

.67D-07 

2 

14D-07 

67D-07 


140-07 

. 66D-07 

3 

. 14D-C7 

66D-07 


. 1 ID- 07 

510-07 

4 - 

. 1 ID- 07 

51D-07 


.34D-15 

440-15 

5 

.340-15 

440-15 


-26D-15 

■310-15 

6 

. 2 63- IS 

• 3 ID- 1 5 


73D-15 

.210-14 

i 

.*>30-15 

.210-14 


43D-15 

. 12D- 14 

e 

430-15 

12D-14 


.4^0-16 

.180-16 

9 

■47D-16 

. 18D- 1 6 


22D-18 

.810-18 

10 

22D- 18 

. 810-18 


■>50-20 

.380-19 

11 

• 7SD-20 

•38D-19 


320-21 

■ 66D-21 

12 

32D-21 

. 66D-21 


■ 160-21 

.550-21 

13 

. 16D-21 

■ 55D-21 


• 48D-22 

.170-21 

14 

460-22 

.170-21 


■360-25 

.580-24 

15 

36D-25 

.580-24 


. 36D-25 

.100-24 

16 

■ 360-25 

.100-24 


. 1BD-25 

. 200-24 


. 29D-25 

.280-24 


• 220-25 

.270-24 


. 3 ED- 2 5 

.260-24 


. 2 8D- 25 

.500-25 

17 

. 280-25 

.500-25 


.360-25 

.700-24 

18 

. 36D-25 

•7CD-24 


28D-25 

. 95D-25 

19 

. 26D-25 

.950-25 


. UD-07 

.760-08 


. 1 1D-Q7 

. 76D-0B 

.000 

.670-07 

.140-07 


67D-07 

.140-07 

.001 

.670-07 

.140-07 


.660-07 

.13D-07 

.008 

660-07 

. 13D-07 


. 510-07 

100-07 

.121 

S1D-C^ 

.100-07 


720-12 

.430-10 

1 OCC 

.■’30-12 

. 4<D- 10 


.280-14 

. 58D-12 

1 OOC 

330-14 

430-12 


32D-14 

.77D-11 

o 

o 

o 

210-14 

730-13 


120- 14 

.420-13 

.236 

13D-14 

.500-13 


. 24D-16 

.760-14 

1 . OOC 

.190-16 

.120-14 


.110-17 

. 1 4D-15 

1 . OOC 

.11D-17 

.160-15 


.140-18 

.25D-16 

1.000 

.830-19 

.150-16 


.550-20 

.87D-17 

1.000 

.350-20 

. 33D-17 


. 56D-20 

.610-17 

1 . 000 

.790-20 

.350-16 


. 26D-19 

.100-16 

.449 

.260-19 

.100-16 


.500-19 

.220-17 

1.000 

.51D-19 

.230-17 


.660-21 

.18D-17 

1 .000 

66D-21 

. 980- IS 


.360-19 

.190-1? 

1 . OCC 

.400-19 

. 560-17 

, 4 03 

. 35D-19 

.520-17 

cec 

. 2 3D- 19 

.110-17 

-0C7 

.930-20 

.680-18 

.001 

.580-20 

.170-17 


. 18D-19 

.960-18 

. 0C1 

.170-19 

.950-18 


.960-21 

.330-18 

.103 

.270-20 

. 16D-1 8 



0 

. 3eD*08 

930*02 

1 

3 60*06 

. 890*02 

2 

■36D*06 

890*02 

3 

360*08 

.890-02 

4 

. 360*08 

890*02 

5 

.360*08 

■ 890*02 

6 

.200*09 

890*02 

0 

.370*08 

.890*02 

1 

47D*09 

890*02 

2 

.220*09 

. e9o*c: 

3 

.200*09 

.890*02 

4 

■47D*C8 

.690*02 

5 

450*De 

. 890*02 

0 

.370*08 

69D*02 

0 

.370*08 

890*02 

c 

. 370*08 

. £90*02 

1 

420*08 

890-02 

0 

.37D*CB 

. eSD-C2 

0 

•37D*08 

.690*:: 

1 

.750*07 

890-C2 


SOLUTION OBTAINED AFTER 20 ITERATION STEPS 
SOLUTION DATA 

SWITCHING POINTS 


NA*<E OF FILE FOR SOLUTION DATA - >BND3D . SOL 


It eventually computes the solution to its own criterion, however, it is 
clear that BND3D has verified the MPMM3D solution. 

The information provided by BND3D.EXTRA is arguable 
essential. This file contains data for the switching function and 
Hamiltonian as functions of time. The plot below is a graphical 
representation of what BND3D.EXTRA provides 
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Useful Information 
in BND3D. REINT 


VII .3. Run BND3D with 
Homotopy 



The Hamiltonian is almost zero, and very close to the tolerance 
The jumps in the Hamiltonian at the switching points is a common 
numerical phenomenon. Also very important, note that this 
verifies the assumed switching structure thrust on at the 
beginning, precisely ten switching points, and thrust on at the end 
Finally, note the hump between the fourth and fifth burns, noting the 
location of such humps is often useful in deciding the location of an 
additional burn 


The file “BND3D. REINT” also supplies useful data in the form of a 
detailed trajectory. The complete state and costate is included The 
plot below, a projection of the trajectory' onto the x-y plane, was 
created using the raw data in the “BND3D.REINT” file. 



This tutorial begins with the solution file from the “Convert 
MPMM3D File to BND3D File, Run BND3D;“ tutorial. 

Suppose we try and accomplish this change in one step, by altering 
the “BND3D.GUESS” file. The script (“BND3D. SCRIPT”) is, 

simply: 
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BND3D. GUESS 
0 
1 
0 
0 

ld-4 
Id- 1 0 
100 

BND3D. SOL 
1 

BND3D. REINT 

Here is the BND3D output to the screen: 


E C S 0 

FREE FINAL TIME : 1 
HOMCTOPY 0 

MU* 1 - OOOOOOOOOOOOOCOOO 

REO* 0 . OCOCCDOCCOCOCOOCCOE*00 

J2* 0 . OOCGOOCODOQOOQODOOE*00 

GO* 1.00000000000000000 

BETA* o OOOCOCOOOCOOGOOOQOE*OC 

R0* 0. OCOOOOOOOOOOOOOOOOE-OO 

RDU* C . COCOCCOOOCGCCCOOOOE*OC 

s* 0 . OOOOCOOOCOCOOODCOQE *00 

CD* 0. OOOOOOOOOCOOOOOOOOE+OO 

ISP* 0 . 567309999999998962 

THRUST* 0 . 516563000000001014 

AI* 1.00000000010536792 

ei* o . ooooooocooooooooooe-»oo 

OMEGA I ■ 0 . OOOOOOOOOOOOOOOOOOE*00 

RA:* 09 . 99999999970668B0 

1-1* 26.5000000009010619 

AF* 6-60014999841091043 

EF* 0 . OOOCOOOOOOOCOOOCOCE+OO 

OMEGAF* 0 . 00000 0000000000 OOOE*00 

RAF* 0 . OCOCOOOOOOOOOCCCOOE*00 

i-F* o . oooocoooooooooocooe*oo 

NOTE; ANGLES MUST BE IN DEGREES 
M* 44 

•N* 25 


INITIAL DATA 


N=25 M*4 4 MS-10 

PRESCRIBED RELAT IVE PRECISION .10D-09 

MAXIMUM PERMITTED NUMBER OF ITERATIONS! 0 0 


IT 

0 

1 

2 

3 

4 

5 

6 


ABS.ERR. 

LEVEL 1 

LEVEL2 

. 15D-02 

. 15D-02 

- 2 2D- 01 

. 15D-02 

. 15D- 02 

. 22D-01 

. 15D-02 

-15D-02 

. 19D-01 

. 15D-02 

. 15D-02 

■ 19D-01 

. 15D-02 

. 15D-02 

. 46D-01 

. 16D-02 

. 1BD-02 

. 22D*00 

. 15D-02 

-15D-02 

. 49D-01 

. 15D-02 

. 15D-02 

. 37D-00 

. 15D-Q2 

. 15D-02 

. 41D*00 

. 15D-02 

. 15D-02 

.90D-01 

-15D-02 

. 15D-02 

.8BD*01 

. 1 5D-02 

.15D-02 

. 13D*02 

■21D-02 

. 2 6D-02 

. 1 8D*02 

. HD-02 

. 1SD-C2 

. 13D-C2 

. 14D-02 

. 15D-02 

. 18D*02 


LEVEL3 

RELAX 

NEW 

CON’D (M) 

. 11D* 03 
. 11D+03 

.000 

0 

. 39D* 08 

. 1 0D*03 
. 1 0D*03 

.003 

0 

. 3*'D*C8 

. 11D*03 
. 13D*03 
. 1 ID* 03 

.026 

.005 

1 

. €7D*06 

. 1 ID* 03 
. 11D*03 

.017 

0 

.370*06 

. 46D*03 
. 4 5D*03 

.011 

0 

. 36D*06 

. 62D*03 
. B2D-03 
. 6 1D*03 

.056 

.013 

0 

. 37D*08 

B1D*03 


0 

.37D*06 


NORM/'M) 
. 93D-C2 

89D*C2 

. B9D-* 02 

. 69D*02 
. 69D*02 
. 69D-02 

89D-32 
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■26D-02 

■36D-02 

.280*02 

.120*04 

C“3 





.140-02 

. 1 5D-C2 

.170*02 

.790- 03 

015 




1 

. 140-02 

. 15D-C2 

.230*02 

. 1 OD* 04 


c 

3 70-* Z 6 

c 


- 3 ID- 02 

.44D-C2 

.400-02 

.170*04 

052 





.140-02 

. 14D-02 

.22D-02 

.100-04 

.019 




6 

■ 1 4T- 0* 

.15D-C2 

.280*02 

.130-04 


0 

360- 06 

c c ~ - 


670- \,2 

.150-01 

.900-02 

.400-04 

.144 



* - - • * a 


. 130-02 

15D- 02 

•27D-C2 

.120*04 

■ 023 




9 

. 1 3D- 02 

.150-02 

.360-02 

.160*04 


o 

.360*06 



5CD-02 

.110-01 

840-02 

.3 60-04 

.155 





13D-02 

.15D-C2 

.34D-02 

1 5D* 04 

.031 




10 

13D-C2 

. 15D-02 

.460-02 

.2 ID* 04 


0 

. 4 00'* O' 6 

Ct* * “ 


• 140-C: 

.300-01 

.220-01 

.700*04 

.135 



C 7*"* w 4. 

11 

. 14D-C1 

.30D-0: 

. 300- C3 

.140*05 


o 

■ 290*0 6 



. 34D-01 

•76D-C: 

.610-03 

.270*05 

.036 



& - 


. 14D-01 

.300-01 

.300-03 

.140*05 

. 004 




(mar-y 

lines orruttea 

fcr brev 

ity] 






62 

-150*00 

. 150*01 

.290*01 

. 24D* 03 


0 

• £20-06 



. 140*00 

.150*01 

.260*01 

.200*03 

089 



c t - 2 

63 

. 140*00 

.i5D*o: 

.260-01 

.170*03 


2 

. 560*06 

c t " - 


. 52D- 01 

. 2 6D*02 

.100-04 

.230*06 

1 .000 





- 140*00' 

.140*0: 

.240-01 

• 150-03 

.062 




64 

. 14D-00 

.140*01 

.230-01 

.260-03 


c 

■ 660*06 



.160-00 

.160-01 

.230*01 

.190*03 

.135 



9zZ'~Z2 

65 

. 1 6D-0C 

.160*01 

.250*01 

.230*03 


o 

730-06 



.110*02 

.300*02 

. 1 CD* 03 

.500-05 

: . coo 



9iZ ~ C i 


•150-00 

.150*01 

.220-01 

. ieo-03 

. 1 01 




66 

.150*00 

150*01 

.240-01 

. 16D-03 


1 

• 920* 06 



. 960*01 

25D*02 

.470*03 

. 560-05 

1 . DOC 



y ^ * 


. 170* OC 

150-01 

.220-01 

.140-03 

.091 





Execution was terminated early because BND3D was clearly stuck 
In this type of situation, where BND3D has difficulty, it is often 
useful to resort to homotopy. 


BND3D has a homotopy loop and is utilized, for this tutorial with the 
following script (supplied as “BND3D 
HOMOTOPY/BND3D. SCRIPT")- 


BND3D. GUESS 
0 
1 
1 

10 

0D0 

0. 516563D0 

0.5673D0 

6.6D0 

0D0 

0D0 

1D0 

ODD 

0D0 

28.5D0 

0 

ld-4 
Id- 7 
100 

BND3D.SOL 

1 

BND3D . REINT 

To make convergence easier, the tolerance was reduced to 10 7 . Ten 
omotopy steps have been suggested and the final semimajor axis 1 * 
requested to be 6.6. 


The output to the screen is very long for a homotopy run, and is 
the tutoria1 ’ however, it may be found in the file 
B.ND3D HOMOTOPY/screen output.” One the other hand, the 
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"BND3D HOMOTOPY/BND3D. REPORT” file 
homotopy progressed: 


indicates how 


the 


ij, 

0 . 

1 . 

■ 2 , 

3 , 

4 , 

5 , 

6 , 

7 , 

8 , 
9 , 

10 f 
11 , 
12 , 
13 . 


25 , 
22 . 

25 , 
30 , 
19 , 
19 . 

26 . 
19 , 
28 . 
-4 , 
17 , 
14 , 
23 , 
17 , 


U. DU 

. 9000000D+00 , 
•8000000D*00 
.7000000D+00 
.60000000*00 , 
•5000000D*00 , 
. 4000000D+00 , 

. 30000000+00 
. 2000000D*00 
. 1000000D*00 , 

. 138777 9D- 1 5 , 
.75000000-01 , 

. 5000000D-01 , 
.25000000-01 , 

. 14 57 16 8D- 1 5 . 


- . 100000 OD* 00 
- . 1000000D* 00 
- . 10000 0 OD*00 
- . 1 0 000 0OD*OO 
- . 10000000*00 
- . 100GDOCD*00 
“ . 1000000D*00 
- . 1 OOQOOOD* 00 
- . 1 0000000*00 
- . 1000000D*00 
- .2500000D-01 
- . 2500000D-01 
- .2500000D-01 
- . 25000 OOD- 01 


s indicates that even though ten steps were suggest, thirteen were 
required. Iterations failed for the ninth step. BND3D then adjusted 
the step size (DU) to one-quarter and continued until completion 


Vn.4. Using MBCM3D 


The following sample input file has been supplied for MBCM3D 

( Tutorials/MBCM3D/MBCM3D. GUESS”): 


Ku* 

Feq= 

J2 = 

Go* 

Beta* 

Rc* 

Rou* 

S* 

Cd= 

isp* 

Thrust* , 

Ai * 

• 1 = 

omega: = , 

RAi * 
i-i* 

Af« 
ef * 

omega f = , 

RAf = 

1-f* 

X ■ 

y * 

z * 
u • 

v * 

W * 

M * 

Lam-X* , 
Lam-Y* , 
Lam-Z* , 
Lam-U* , 
Lam-V* , 
Lam- W* , 
Lam-M* , 

TF r 
G1 ■ 

G2 * 

G3 * 

G4 - 
G5 * 

G6 = 

G7 < 

G8 * 

G9 . 

G10* 

ACC* 

ITM= , 100 


: -o::::cococccocs 

c . OOOCCOOSOOOOOOO 
c . occcoc-ccccooooo 

0 OC960COGCOOOOOQ 
o. occoc::cooo::oo 
0 . COOOGOOOCDOOCCO 

c.ooococoooooocoo 
o . occo:::o:oocc:o 
c . ocooccoccooccoo 

134 . OOCCOOOCCOOCOCO 
0 0300CCDOOOOOOOO 

3 8 4 7 3 05 000000000 
0 . 023777042000000 

0 oooocooooococoo 
o . ocooocoocoooooo 

0 . 000C000000C0C 00 
1 ■ 5000C00000CC000 
0.333333333333333 
0 . cocccocoooococc 
c. oooooooocoooocc 
c . occoooccooooooo 

*3 . 1376eC190F73i56 
2 ■ 3”50C7893528269 
o . oocoooooooooooo 
-c .309133504323169 
-0.393443660534349 
o . ocooocoocoooooo 
1 . 527GOCOGOOQCOCO 
0 084150649480784 
-C.C70063915C70165 
0 . 000000000000000 
0.531758699754281 
0.737783173534899 

c . ooooooooocooooo 

0 782111317020586 
19 . 053149861397220 
0 00000000000C000 
0. 000000000000000 
-0.656087-»95635957 
-0 . 235988651457670 
-0.000453092937198 
0 OOOOGCOCOOOOOOO 
0 000000000000000 
0.205432772910901 
-0 . 0286054*0141037 
0 006351966699377 
0 000001001000C-00 
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note all angles ere in degrees 


The MBCM3D iterations, output to the screen (see file 
“Tutorials/MBCM3D/screen. output") follow: 


ITERATIONS * 1 CALLS OF VF02AD * 3 


X 


-0 31 1*7680 19 0673 2£+01 
-0. 309133504322] 7E*00 
0. 1527COOOOOOODOE-D1 

o . oooooccrooocooE*oo 

0 . OOOOOCOCOOOOOOE-OO 
0. ODCOOCCOOOOCOOE-OC 

-0.2359B865145767E*00 
0 . OOOOOOOOOOOOOOE* 00 
0. 6351966699377 0E- 02 


0 . 23750078935283E* 01 
-0 . 39344366053 43 5E* 00 
0.842 506494 6 07 84 E- 01 
0 . 53175B6997542BE«C0 
0.7B211131702C59E-0C 
o. oo ooooooo oco:oe*c: 
-0. 453092937198001-03 
0.2 05 43277291 090E*0C 


F * 0.86334 4 94474323E-G7 


C 


0. OOOCOOOOOOOOOOE-OO 
0 .3180397B123701E-02 
C.4375169:e96031E*C2 
C . 84402672 955175E- 03 
0 OOOOCOOOOOOOOOE*C0 
0.22530999940806E-03 
-C . 2274086 974 996 0E- 04 
-C. 575695 3 1495894 E- 05 


0.0000 0000000 00 0E- CO 
0 . 1 67 6456 0733 063 E- 02 
0. OOOOOOOCOOOOCOE-OC 
0 . 00000000000 OOOE* CO 
0 . 0000 C 000000 DOCE-OC 
C. 1733 315396964 5E -03 
■0.2;i7739363C278E-04 
0. 9062 3565540662 E-C 5 


ITERATIONS * 2 CALLS OF VFC2AL = 2 


X = -C . 31C365C2209544E-C1 
-C . 31 13636 342 5003E* OC 
C. 15270COCCCOOOOE*C1 
-C 801 980941 0634 0E- 20 
-0.464S1533630C1 5E-20 
-0. 3C:28:7ie97062E-19 
-0.236206 597 57290E*OC 
0.24944301628696E-20 
0. 6465666S122166E-C2 


0.2392B642I4496CE-C1 
-0.3917 53 2 1262 767E-0C 
0.8376 53 141956B0E -01 
0 . 535965621 5523 6E-C0 
0 . 7 82 03 61 8243 91 7E*C 0 
0.12232 03 9786695E- 19 
-0 . 2056 96C5044 466E- 03 
0 . 2055405331371 9E* 00 


f * c. 249 :c:o 7 ieo 5 iiE-i: 


0 . 00001 0 00 DC 00 00*-" 

o.oooc: :ococ:c:ci- ;: 
-0.7C0€35:507C165E-C: 
0.73776 3 173534SCE*c r 
0 . 19052 1498613 57E- C2 
-c 65e:e'7?it 2 ^ 

0 . OCCOOCOCCCOCOOE-CC 
-0.2860 5410141C37E-C1 


-0.226:565232C329E-C2 
-c . 2926 :?“: 909 " 9 :e-:: 
-C -liI1569e66C4! A E-01 
-0.16C25476632695E-:: 
-0. 7426722973491 5E- s 

o . ooc:oc:occ::c:f- 
o . cccc:::occ:::cf- 
o . occoC'C:cc;::c:e-; ■ 


0 . 25568 446552:98*- 15 
-0 4252336262341SE-2: 
-0.70567462291251E-C: 
0.734455690936E6E-CC 
0 . 19053764 9149"'5E* 02 
-0 6566£256:i:73SE.c: 
-C .374;5954:5336’ 7 E-2C 
-0.264156726C4799E-C1 


C * -C . 58042563465467E-2C -C . 30064192465134E-20 
-0 . 868 632300431 7 8E-03 0.627799253 62 639E-03 

-C 7756652 6491 424E-C3 C.66247067359436E-22 

-0. 65654CE22S1S2CE-03 -C . 1 54936 9 £704 C00E- 1 9 
-0 . 156742 19322 67 6E- 21 -C . 211 5 93 653 0661 5E- 19 
-0 . 3 6 5C51 0906 133 EE -04 0 .299112571 19 057 E- 04 

0.362C67 677136 02E-06 -0.8532 57 21599261 E-C 7 
0. 242366 7 06 0626 3E- 05 0. 3067486376860 6E-05 


(lines omitted] 


ITERATIONS * 5 

X = -C. 3111316 62 53579E-01 
-0 .310097 034 1424 0E*0G 
0.1527CCOOODOOOCE-C1 
-C. 55623767494574 E- 20 
-C.79958433660922E-2G 
-C .75769120625379E-20 
-0 . 23 59665475 9671E* 00 
-C. 7977133422306 4E -17 
0 62519253321692E-C2 


CALLS OF VFC2AL = 5 

C. 236306940 92 665E-C1 
-0.392731270 031 03 E*0C 
C. 839490 69 044446E-01 
0.5234663::E974:e-C: 
0.7 62 11132986 574 E-C! 
-0. 156952 0635S46BE-20 
-0.45317463C73C11E-03 
C.20543273853632E-0C 


F * 0.212172 02 906824E-25 


C . 112766b 6356£“ , 6E-:2 
0 .50606043C77460E-03 
C .1366C752767354E- Cl 
0 . 49671708:'716S2E-03 
C 63852 3E931E222E- :< 

c.coc::c:::o::o:e*:: 

0 . 10191607C7i;iOE-19 
C . 69:50*26954e5SE-20 


o.€95ob:o9E"712ce-:: 
-C. 2 C 264 S“E 2 " 6292 E- 2 C 
-c 703::5s:"3?55£r- :: 

0 .73634 C49255129E-C0 
0 :?072766:66345E-:: 
-0 . 6560S762136251E*:: 
-0 .626}e395;C65£fr-:- 
-0.26605461B566I7E-Ci 


C 


- 0.5812999 87 686B2E- 20 
-0.77946793908268E-07 
-0 . 711327 6B54 S312E- 07 
-0.496519123 59394E- 07 
-0.20B84662B8415BE-20 
-C.49293613635371E-09 
C. 5140332604014 5E -11 
0.33 07221163 5954E- 10 


0. 89436352481144E-21 
0 . 3301 8123 6851 69E-06 
0. 6725466625966 9E -20 
0 . 4 60734029654 87E- 20 
-0.84790 04 31647S7E- 20 
0.35773 77 2632925E-0S 
-0. 17852386235973E-11 
C. 2116706 6096293E-10 


0.20736C66980633E-06 
0.20C01798439‘ t 50E-C^ 
-0. 3590747965621 6E-C7 
0.109 04 6612 64 651E-06 
0.8312E526264761E-09 
0.71C54273576C: 01 -14 
-0. 34634 0337C 02 15E-21 
0 . 396023255 :3391E- 16 


THE PRINTING OF THE LAST ITERATION GIVES THE 
VALUES THAT ARE RETURNED BY SUBROUTINE VFC2AD 


SOLUTION CONVERGED--- 

X * -0.311131€62535' > 9£*D1 

0.2363 08 9409266 5E*C1 
0 . 699 08 1 09877 12 0E -20 
-0.310 09 70341424 OE *00 
-0.39273127003103E-00 
-C.2C264578276292E-20 
0.15270000000000E*:i 
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0 83949089044448E-01 
-0. 703 0059 0"39595E-0J 
-0 5 582 37 674 94 57 4E- 20 
0. 53348632269741 E* 00 
0.73634049255129E*00 
-0.79958433860932E-2D 
C.76:ii: 32986S74E«0C 
0 . 19C73788066349E-02 
-0 .75789120B25379E-20 
>0 156 952063 59 4 68E- 20 
>0 656 08762138191 E* 00 
-0.23598854759671E-00 
-0 45317463073CilE-03 
-0 62836395106568E-n 
-C 79771334223064E-17 
0.20543273853632E*00 
-0.286D5461858627E-01 
0.63519253321692E-02 
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Appendix A GSHOOTs File Format 

The input file - for “GSHOOT” has a specific file format The file 
must consist of exactly 14 lines. The variables read from this file 
have a specific order: MU, GO, ISP, THRUST, MO, AO, EO, WO. 
AD, ED, WD, TMAX, NGS, and NDO All variables are of the type 
REAL except the last two, NGS and NIX, which are of the type 
INTEGER An example file is listed below. 

Ku * ] CC 

&C * 3 .CO 

3sp a 0.56^3 

Thrust * 0.5166 
He c iz CCOC 
• ~ * 1.00000 
eo r C.C00 
V? * C.00C 

* 1 2P5 

* 0.219 
wc * 0. 000 

tkax = c.occ 

NOS * IOC 
NIX * 3 

On each line intended to supply a REAL variable, the FORTRAN 
FORMAT layout is (1X,A9,F30.15); for INTEGER variables, this 
statement is (1X,A9,II0). Therefore, each line starts with a blank 
space followed by nine characters, all of which are ignored. Only 
the numerical data following is used. 
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Appendix B The PAT2D and PAT3D File Formats 

The PAT2D file format is used by MPMM2D and PAT2D. The 
PAT3D file format is only used by MPMM3D They are called the 
PAT formats because all of the information supplied by the PAT2D 
format is used by PAT2D; only some of the information is used by 
MPMM2D and MPMM3D. Exactly what information is used by 
MPMM2D and MPMM3D is described in Chapter IV. 


The PAT2D format is represented below: 


TCL * I 

MU ■ • 

T * A 
GO - • 
ISP « • 

AO * I 
EXO • • 
EYO * « 
PS > • 

DOF * • 

nT . • 
a t « 
ex * • 
ey * • 

NORB « 2 

NODE * 3 

SEL * 1 


ndex , x 

y . u . v 

rr , 

lx. ly, 

lu, Iv, 

lift 

tf , 

91.92 . 

93 . 

94 

9 - , 06 

1 . • , 

• , • , 

$ , 

9 . «, 

#, * . 

• 

• . 


1 . 

• 

• | 

2, • , 

«, *, 

« . 

*, #, 

1, #, 

• 

1 . 

# «, 

1, 

• 

• . f ’ 

3 , • 

i • , 

« . 

* • . 

* . * . 

t 

# 

# 1, 

• . 

• 

9 , 9 . 

4 . • 

* * 

• « , 

« , 

* . * . 

• 1, 

« . 

• . 

• • , 

• 

« 

9. #, 


ex = a 


ey * « 

NODE = 3 

SEL = 2 


1 r.aex . x 

y ' 

u . v , in , 

lx , ly, 

lu. 

lv, lift, tf , gl . 

92,93, 

94 g6 

i . * 

• , 

• , * . 

• , 9 . 

« 

» », i , 1 , 


A, A A « 

2. 9 . 

• 

9 . «, 

9 . 9 

• 

«. * « «, 


I, * 9. 9, 

1 , 

• . 

i . 

* « , 

9 9 . 

* . 

#, 9 , | « 

* , * , 

9. • , •. 9. 

9. 9. 9, A. 

4 

A 

i . 

* « 

• , 

1 • . 

9. 1, 

A , 

• •, •, *, 

9. 9, 


s f 
* « 


ex 

NODE = 3 

SEL * 3 

INDEX.X, Y,U.V.K.TF,L1.L2 


1. 

2 . #, 
3 . •. 

4 , $ , 


i . « 


where the symbol "#" is used in place of digits. The first eleven 
lines give constants for the orbit transfer problem in type REAL 
These have a fixed order. TOL, MU, T, GO, ISP, AO, EXO EYO AF 
EXF, and EYF. Their descriptions follow: 


TOL . 
MU . . 
T . . . 
GO . . 

ISP . 
AO . . 
EXO . 
EYO . 
AF . . 
EXF . 
EYF . 


THE SOLUTION TOLERANCE 

THE GRAVITATIONAL CONSTANT FOR THE CENTRAL BODY 

THE THRUST LEVEL OF THE ROCKET MOTOR 

EARTH'S GRAVITATIONAL ACCELERATION AT SEA-LEVEL 

[ONLY USED FOR GET MOTOR FUEL CONSUMPTION] 

SPECIFIC IMPULSE OF ROCKET MOTOR 

INITIAL ORBIT SEMIMAJOR AXIS 

INITIAL ORBIT X-COMPONENT ECCENTRICITY 

INITIAL ORBIT Y-COMPONENT ECCENTRICITY 

FINAL ORBIT SEMIMAJOR AXIS 

FINAL ORBIT X-COMPONENT ECCENTRICITY 

FINAL ORBIT Y-COMPONENT ECCENTRICITY 
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Note that these apply to the transfer as a whole, esp. when referring 
to the initial and final orbits. The FORTRAN FORMAT edit 
descriptors for each of these first eleven lines is (1X,A6,E27.20) 

The PAT3D format up to this point is identical except that HXO, 
HYO, HZO, EXO, EYO, HXF, HYF, HZF, EXF, EYF replace AO, 
EXO, EYO, AF, EXF, and EYF. Their descriptions follow: 


HXO INITIAL ORBIT X-COMPONENT ANG . MOMENTUM 

HYO INITIAL ORBIT Y-COMPONENT ANG. MOMENTUM 

H20 INITIAL ORBIT 2-COMPONENT ANG. MOMENTUM 

EXO INITIAL ORBIT X-COMPONENT ECCENTRICITY 

EYO INITIAL ORBIT Y-COMPONENT ECCENTRICITY 

HXF FINAL ORBIT X-COMPONENT ANG. MOMENTUM 

HYF FINAL ORBIT Y-COMPONENT ANG. MOMENTUM 

H2F FINAL ORBIT 2-COMPONENT ANG. MOMENTUM 

EXF FINAL ORBIT X-COMPONENT ECCENTRICITY 

EYF FINAL ORBIT Y-COMPONENT ECCENTRICITY 


For both PAT2D and PAT3D formats, the next line indicates how 
many intermediate transfer orbits there are. The variable NORR 
takes on this value. The FORTRAN FORMAT edit descriptors for 
this line is (1X^6, 13). This same layout is used for the next two 
lines, both also containing INTEGER data. These lines specify 
data for the first bum. NODE is how many nodes, not counting the 
first one, are to be used for this bum. Specifying a “3” for NODE 
indicates that four lines of data will describe the burn. 

The line after NODE’S is for SEL. The variable SEL indicates 
which method should be used. Note that in the PAT2D 
representation above, three different values are given for SEL. A 
“I” indicates that the data below is in a multiple-point shooting 
format but Direct Collocation with Nonlinear Programming 
(DCNLP) should be used in the first attempt to obtain a solution. A 
“2” also indicates that the data below is in a multiple-point shooting 
format but that multiple-point shooting should be used in the first 
attempt to obtain a solution A “3” indicates that the data below is in a 
DCNLP format and DCNLP should be used in the first attempt to 
obtain a solution. The following table summarizes: 


SEL 

Guess Format 

Method to try First 

1 

Multiple Shooting 

DCNLP 

2 

Multiple Shooting 

Multiple Shooting 

3 

DCNLP 

DCNLP 


No matter what format the data lines will be in, the line following 
SEL’s line has the FORMAT edit descriptors (IX, A). The contents of 
this line are ignored. 

Note that since MPMM3D cannot accept SEL=3, in PAT2D only 
SEL=1 or SEL=2 is acceptable. 

The next NODE+1 lines are the guess data for that bum. The 
FORMAT edit descriptors are (IX, 13^1,50(027. 20^1) irrespective 
of which guess format is intended. Considering only PAT2D, the 
multiple-point shooting format has 18 elements in each line. These 
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elements wcu, the following order: INDEX, X, Y, U, V, M, LX LY 

Hne MLhp fi 'J * Gl * G2 ’ G3, ° 4, G5 ’ 06 “ INDEX ” numbers each ’ 
line, the first line represents the initial point for this burn and last 

line represents the final point for this burn. The lines for each burn 
are evenly spaced. % Y, U, V are the Cartesian components of the 
2D posiuonandvelocity vectors, respectively. “M"is?hemass 

,“*’, LY ’ LU ’ LV ’ LM are the values of the Lagrange multiplier 
functions/ or costates, V X v , and A m , respectively. “TF” is the 

length of time the burn lasts. “Gl, G2, G3” are the constant 
grange multipliers, v f , associated with the final boundary 
conditions. “G4, G5, G6” are the constant Lagrange multipliers v 
associated with the initial boundary conditions. ’ °’ 

ra°cMfn? D ^eL m p^ tiPle ; P ° int Sh °f ing format has 26 "Aments in 
7 IT v w arC m the followin g o^er: INDEX, X Y 

« b:' G8 G^'lO 'if' LU ’ LV ' LW ' LM - TF ' G1 ' M. 04 cl 

from PAT21? ’ Their meanings are simple extensions of those 

pe DCNLP format has 9 elements in each line. These elements 
are in the following order: INDEX, X, Y, U, V M TF Ll L2 All of 
these are as described above, except “Ll,’ L2" w^h L the Cart^iaL 
mponents in the inertial frame of the thrust direction unit vector 
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